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SUMMARY 


OBJECTIVE 

Cox  and  Munk  (1954)  proposed  that  the  fluctuating  surface  of  the  open  ocean  could  be 
accurately  represented  by  many  small  mirror-like  facets,  each  of  which  is  randomly  tilted  with 
respect  to  the  local  horizon.  Through  aerial  photographs  of  sun  glint,  they  determined  the 
statistical  distribution  of  capillary  wave  slopes  as  a  function  of  wind  velocity.  However,  when 
their  equation  connecting  the  slope  distribution  with  sun  glint  is  used  on  the  horizon,  an  infinite 
glint  is  predicted  even  though  it  can  easily  be  shown  that  sun  glint  should  never  exceed  solar 
radiance. 

The  objective  of  this  report  is  to  remove  the  infinite  horizon  radiance  from  the  Cox-Munk 
model  of  the  ocean  surface. 

RESULTS 

An  integral  equation  for  ocean  radiance  is  derived  in  this  report.  This  equation  predicts  a 
finite  radiance  on  the  ocean  horizon  and  the  proper  value  of  sun  glint  in  the  calm  sea  limit.  Away 
from  the  horizon,  an  approximation  can  be  used  that  reduces  this  equation  to  the  form  used  by 
Cox  and  Munk.  The  unphysical  infinity  resulted  from  improper  use  of  an  approximation  in  a 
region  where  it  is  not  valid. 
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1.  INTRODUCTION 


Forty  years  ago,  Cox  and  Munk  (1954, 1956)  presented  a  model  for  the  fluctuating  surface 
of  the  open  ocean.  They  supposed  that  the  sea  surface  consisted  of  numerous  wave  facets.  Each 
facet  is  a  flat  surface  with  an  area  perhaps  as  small  as  a  few  mm^.  At  any  given  instant,  the 
ocean  surface  can  be  regarded  as  a  collection  of  facets  that  are  randomly  tilted  with  respect  to 
the  local  horizon.  As  time  passes,  the  tilt  of  a  given  facet  varies  under  the  influence  of  the  wind. 
When  the  solar  disk  is  reflected  by  the  open  ocean,  these  fluctuating  facets  produce  a  dancing 
pattern  of  highlights  known  as  sun  glint. 

By  photographing  solar  reflections  from  an  airplane,  Cox  and  Munk  were  able  to  measure 
the  probability 

P  =  p(s^,Sy,w)dS^dSy,  (1) 

that  a  wave  facet  would  have  a  slope  within  ±  dSx/2  of  Sx  and  ±  dSy/2  of  Sy.  Sx  and  Sy  are  the 
slopes  of  the  facet  in  the  upwind  and  crosswind  directions  respectively,  and  W  is  the  windspeed. 
Cox  and  Munk  obtained  the  following  expression  for  the  wave  slope  probability  density: 


^  C?  ^  ^  PVD 1 

» 

1 

f  7 

1 

"*2 

1  “ 

^  A 

CT2=0.000  +  3.16  10"^r 
u 


CT^  =  0.003  + 1.92- 10“^  fF. 
c 


(2) 


This  is  actually  the  lowest  order  term  of  the  fit  found  by  Cox  and  Munk,  but  it  will  be  accu¬ 
rate  enough  for  our  purposes.  Here  and  Oc  are  the  standard  deviations  in  the  upwind  and 
crosswind  directions,  respectively.  Both  depend  on  windspeed,^  which  must  be  entered  into  these 
expressions  in  m  s~^. 


The  Cox-Munk  model  is  often  cited  in  oceanographic  literature  and  has  proven  useful  in  the 
forty  years  since  it  was  first  introduced.  An  interesting  history  of  the  Cox-Munk  contribution  and 
its  relationship  to  earlier  measurements  made  by  Duntley  was  given  by  Preisendorfer  and 
Mobley  (1986). 


The  primary  goal  of  the  original  Cox-Munk  work  was  to  determine  the  wave  slope  probabili¬ 
ty  distribution,  previously  mentioned.  Now  that  the  wave  slope  distribution  is  well  accepted,  the 
concern  is  no  longer  with  deriving  the  wave  slope  probability.  Instead,  we  are  interested  in  the 


*  A  windspeed  of  1  m  s~*  is  typical  of  “light  air;  on  land  the  wind  direction  is  shown  by  smoke  drift  but  wind  vanes 
do  not  move  and  at  sea  ripples  with  the  appearance  of  scales  are  formed,  but  without  foam  crests.”  A  windspeed 
of  10  m  s~^  is  typical  of  a  “fresh  breeze:  on  land  small  trees  in  leaf  begin  to  sway  and  at  sea  there  are  moderate 
waves  of  1-Va  to  2-V2  m  height,  many  whitecaps,  and  some  spray.”  (Bowditch,  1984). 
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inverse  situation  where,  given  the  value  of  p,  we  wish  to  predict  the  radiance  reflected  from  the 
surface  of  the  sea.  In  particular,  we  are  interested  in  measuring  the  reflected  solar  radiance  from 
a  low  altitude,  such  as  the  deck  of  a  ship,  when  the  sun  is  anywhere  in  the  sky,  including  setting; 
we  want  to  look  straight  into  a  bright  pattern  of  solar  glint  on  the  horizon  and  predict  its  radiance 
at  a  given  wind  velocity. 

This  prediction  has  already  been  provided  by  Cox  and  Munk  (1954).  Their  equation  (9)  gives 
the  radiance  N  leaving  a  glint  pattern  when  the  solar  irradiance  arriving  at  the  glint  pattern  is  H. 
In  rearranged  form  and  our  notation,  their  prediction  is 

—  =  - sec0;.psec  d„p  (Cox  and  Munk,  1954  (9)). 
li  4 

In  this  equation,  0;.  is  the  zenith  angle  of  the  reflected  ray,  p  is  the  reflectivity  of  sea  water, 

0„  is  the  tilt  of  the  facet,  and  p  is  given  by  equation  (2).  The  geometry  implicit  in  this  equation  is 
as  follows;  the  incident  ray  comes  from  the  center  of  the  sun,  the  reflected  ray  goes  to  the 
receiver,  and  the  facet  slope  is  such  that  a  specular  reflection  occurs  between  the  two  rays. 

A  problem  exists  with  the  equation  (Cox  and  Munk,  1954  (9))  for  the  particular  geometry  in 
which  we  are  interested.  When  viewing  the  horizon,  0^  approaches  90°  and  the  radiance  pre¬ 
dicted  by  Cox  and  Munk  (1954  (9))  approaches  infinity.  But  an  infinite  radiance  is  never 
observed.  A  perfectly  calm  ocean  (which  can  never  occur  but  serves  as  a  useful  limiting  case), 
which  reflects  the  rays  of  the  setting  sun  into  a  receiver  mounted  at  the  edge  of  the  shore  can  be 
imagined.  Then,  any  theory  of  solar  glints  should  predict  a  radiance  equal  to  the  solar  radiance 
times  the  reflectivity  of  sea  water,  because  a  perfectly  calm  ocean  is  exactly  like  a  flat  mirror. 

This  infinity  is  the  problem  addressed  by  this  report.  Three  possible  classes  of  solution  are 
anticipated. 

(1)  The  theory  is  incomplete.  The  effects  that  have  been  left  out  are  responsible  for  the 

infinity. 

(2)  The  theory,  in  particular  the  equation  (Cox  and  Munk,  1954  (9))  above,  is  wrong.  The 

infinity  results  from  the  error. 

(3)  The  theory  is  correct  over  a  restricted  range  of  conditions.  The  infinity  results  from  the 

application  of  the  equation  (Cox  and  Munk,  1954  (9))  outside  this  range. 

In  returning  to  Cox  and  Munk  (1954  (9))  and  trying  to  understand  its  derivation,  we  bene¬ 
fited  from  a  paper  by  Plass,  Kattawar,  and  Guinn  (1975),  who  verbally  restated  the  geometry 
surrounding  Cox  and  Munk  (1954  (9))  in  such  a  way  that  we  could  derive  our  own  version  of 
Cox  and  Munk  (1954  (9))  instead.  This  report  presents  that  derivation.  The  resolution  of  the 
unphysical  infinity  in  Cox  and  Munk  (1954  (9))  was  discovered  to  reside  in  class  (3)  above. 
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2.  THE  GEOMETRY  OF  REFLECTION 

In  this  report,  we  will  only  consider  single  scattering  events  as  they  occur  at  the  pixel  foot¬ 
print.  Multiple  scattering  and  shadowing  within  the  footprint,  and  transmission  losses  and 
refraction  as  the  beams  traverse  the  atmosphere  to  and  from  the  footprint  will  all  be  ignored. 
Polarization  will  be  ignored.  Gravity  waves  will  not  be  included. 

Figure  1  shows  the  geometry  of  reflection  and  introduces  the  symbols^  that  will  be  used  by 
this  report.  We  have  chosen  a  coordinate  system  whose  origin  is  the  point  of  reflection  with  the 
X  axis  pointing  upwind,  the  Z  axis  toward  the  zenith,  and  the  Y  axis  crosswind  such  that  a  right- 
handed  system  is  formed.  The  X-Y  plane  is,  therefore,  horizontal  at  the  point  of  reflection.  The 
tilted  facet  passes  through  the  origin.  Three  unit  vectors  u  «  (0, 4))  are  involved  in  reflection. 
Each  unit  vector  has  polar  coordinates  0,  the  zenith  angle,  and  (J),  the  azimuth.  Two  of  these  unit 
vectors  are  shown  in  figure  1:  U5,  pointing  from  the  origin  to  the  source,  and  iv  pointing  from 
the  origin  to  the  receiver.  Note  that  the  zenith  angle  of  the  facet  normal  (not  shown)  is  the  same 
as  the  tilt^  of  the  facet. 


I 


Figure  1.  The  geometry  of  facet  reflection.  Unit  vectors  Ug  and  Ur  point  toward  the 
source  and  receiver,  respectively,  from  the  origin  of  a  right-handed  coordinate  system 
located  at  the  reflection  point.  The  X  axis  points  upwind,  the  Y  axis  points  crosswind, 
and  the  Z  axis  points  toward  the  zenith.  The  facet  normal  Up  has  been  left  out  of  the  fig¬ 
ure  for  clarity.  Azimuths  are  considered  positive  when  measured  in  the  counterclock¬ 
wise  direction  looking  down  along  the  nadir,  a  convention  that  conforms  with  mathemat¬ 
ical  usage  but  does  not  conform  with  nautical  usage. 


2  The  symbol  list  at  the  end  of  this  report  correlates  our  symbols  with  those  used  by  Cox  and  Munk. 

3  The  tilt  is  the  angle  of  steepest  ascent  within  the  facet. 
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The  three  vectors  for  source,  facet  normal,  and  receiver  are  connected  by  the  law  of 
reflection; 

+  U^  =2c0Sfi)U^  ^  (3) 

where  o)  is  the  angle  of  incidence  and  the  angle  of  reflection. 

3.  THE  OCCURRENCE  PROBABILITY 

We  call  p  the  occurrence  probability  density  because  it  gives  the  likelihood  that  a  facet  with  a 
given  slope  will  occur  in  a  group  of  facets  with  all  possible  slopes.  The  shape  of  p  is  shown  in 
figures  2  and  3  for  windspeeds  of  1  and  10  m  s“^  respectively.  The  volume  of  a  column  under¬ 
neath  p,  whose  base  is  dSx  on  one  side  and  dSy  on  the  other,  is  the  probability  P  given  by 
equation  (1).  The  entire  volume  under  the  surface  of  p,  which  is  the  same  as  the  integral  over  all 
slopes,  is  unity. 


z 


Figure  2.  A  plot  of  the  occurrence  probability  density  pas  a  function  of  slope  for  a 
windspeed  of  1  m  s~^  In  the  calm  sea  limit,  this  function  approaches  a  two-dimensional 
delta  function  selecting  zero  slope  in  each  direction.  The  coordinate  system  of  figure  1 
has  been  inserted  at  the  top  of  the  figure  to  illustrate  the  relationship  between  coordi¬ 
nates  and  slopes  put  forth  in  equation  (A5).  Note  that  the  first  quadrant  contains  nega¬ 
tive  slopes. 

By  using  the  occurrence  probability  density,  a  specified  area  of  the  ocean  surface  can  be  used 
to  represent  the  ocean  as  a  whole.  Let  this  representative  piece  of  the  ocean  have  an  area  A.  Area 
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A  can  also  be  thought  of  as  the  footprint  of  a  single  pixel  in  an  image  of  the  open  ocean.  Within 
the  footprint,  a  typical  facet  with  area  dF  and  tilt  Q„  is  shown  in  figure  4.  Choose  a  convenient 
small  range  of  X  and  Y  slopes,  which  will  be  the  same  for  all  facets.  Let  the  horizontal  projection 
of  each  facet  occupy  an  area  dA  whose  ratio  to  the  total  area  A  is  equal  to  the  probability  that  the 
facet  will  occur: 

dA 

—  ^pdS^dSy.  (4) 

Since  the  occurrence  probability  density  has  been  normalized  to  unity  when  integrated  over 
all  possible  slopes,  we  have  the  sensible  relationship  that  the  footprint  area  is  the  sum  of  individ¬ 
ual  horizontal  facet  projections: 

\dA  =  A\\pdS^dSy=A,  (5) 

The  collection  of  all  facets  similar  to  the  one  shown  in  figure  4  will  represent  the  average 
condition  of  the  wind-driven,  open  ocean  for  any  chosen  wind  ''elocity.  These  facets  can  be 
connected  together  to  make  a  continuous  surface  whose  shape  uas  been  studied  by  Gordon 
(1969)  and  Preisendorfer  (1976),  They  called  this  surface  the  “ergodic  cap.” 
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Figure  3.  A  plot  of  the  occurrence  probability  density  p  throughout  slope  space  for  a 
windspeed  of  10  m 
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Figure  4.  Upper  figure:  A  representative  facet  of  area  dF  from  the  pixel  footprint.  Its  pro¬ 
jection  dA  onto  the  horizon  is  proportional  to  the  probability  of  finding  this  facet  within 
the  footprint  Its  projection  dB  in  the  direction  u  gives  the  area  of  the  facet  seen  from  u 
and,  hence,  gives  the  relative  importance  of  this  facet  with  respect  to  all  other  facets 
(not  shown  here)  that  populate  the  footprint  Area  dB  is  proportional  to  the  probability 
this  facet  has  of  interacting  with  a  beam  toward  u.  Lower  figure:  Coordinate  system  for 
the  upper  figure.  Unit  vector  Up  is  normal  to  the  facet  During  the  integration  of  equa¬ 
tions  (8),  (12),  and  (13),  u  remains  fixed,  Up  ando)  vary,  and  no  sort  of  reflection  is 
considered  to  occur. 
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4.  THE  INTERACTION  PROBABILITY 

The  occurrence  probability  is  the  same  as  the  probability  that  a  vertical  beam  will  interact 
with  the  facet.  However,  as  pointed  out  by  Plass  et  al.  (1975),  a  beam  that  is  not  vertical  will 
interact  with  the  facet  according  to  a  different  probability,  which  we  will  denote  by  Q  and  call 
the  interaction  probability.  The  interaction  probability 


Q  =  q(e,<l,X.Sy.fv)dS^dS, 


(6) 


gives  the  chance  that  a  beam  with  the  arbitrary  direction  u  *  (0,  (J))  will  interact  with  a  facet 
whose  slope  is  within  ±  dSx/2  of  Sx  and  ±  dSy/2  of  Sy  when  the  windspeed  is  W. 

We  will  now  derive  the  mathematical  form  for  q  and  relate  it  to  p.  The  interaction  probability 
between  the  beam  and  the  facet  under  consideration  will  be  proportional  to  the  area  dB  of  that 
facet  projected  normal  to  direction  of  the  beam  as  shown  in  figure  4.  Hence,  by  reasoning 
similar  to  that  used  to  arrive  at  equation  (4),  we  have 

—  ^qdSxdSy  .  (7) 

D 


In  this  equation,  B  is  the  sum  of  all  other  similarly  projected  facets: 


B=  jdB 


u=:const 


(8) 


The  integral  of  equation  (8)  is  restricted  to  those  facets  for  which  the  beam  will  strike  the 
front  rather  than  the  back  of  the  facet;  that  is,  for  all  slopes  such  that  (os  nl2.  Furthermore,  the 
notation  “u  =  const”  has  been  added  here  to  indicate  that  the  beam  direction  is  held  constant 
during  integration  over  the  slopes. 

The  interaction  probability  can  be  expressed  in  terms  of  the  occurrence  probability  because 
of  a  geometrical  relationship  between  dB  and  dA  evident  from  figure  4: 

dB-dFcos(0 


dA  =  dF  cos  Or 


Hence, 


dA 


cos© 


COS0, 


n 


which,  together  with  equations  (4)  and  (7),  implies  that 

A  cos© 

9- - P 

B  cosOr, 


(9) 


(10) 


(11) 


Equation  (11)  is  the  desired  relationship  between  q  and  p.  For  a  beam  directed  to  the  zenith, 
the  angle  of  incidence  is  the  same  as  the  facet  tilt.  Then  (10)  shows  that  dB  equals  dA,  which  is 
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also  obvious  from  figure  4  when  u  points  along  Z.  Therefore,  B  equals  A  and  (11)  shows  that 
q  equals p  when  the  beam  is  vertical. 


Using  equations  (7)  and  (11)  in  (8),  we  arrive  at  the  following  expression  for  the  area 
B(0,  <j),  W)  in  terms  of  the  occurrence  probability: 


5 

A 


JJ 


cos  6? 

COS0„ 


\i=comt. 


pdS^dSy. 


(12) 


B  must  be  evaluated  under  the  exclusion  of  facet  backs  ws  ji/2,  for  each  direction  of  the  beam, 
and  for  each  value  of  the  wind  velocity.  The  behavior  of  BIA  for  a  footprint  of  unit  area  is  shown 
in  figure  5  for  various  wind  conditions. 


Figure  5.  The  area  of  the  ergodic  cap  projected  toward  the  beam  for  various  values  of 
windspeed.  The  beam  is  upwind  (zero  azimuth)  and  a  unit  area  has  been  assumed  for 
the  footprint.  The  footprint  contains  a  surface  that  acts  rough  so  that  it  can  be  seen 
even  from  a  glancing  direction  (zenith  angle  90°).  The  greater  the  wind,  the  rougher  the 
surface.  On  the  other  hand,  the  flat  mirror-like  surface,  which  represents  mean  sea 
level  inside  the  footprint,  has  a  projection  that  varies  as  the  cosine  of  the  beam  zenith 
angle.  It  vanishes  at  a  glancing  angle. 
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Using  (12)  in  (11),  we  have 


03<nl2 


COS  ft) 
COS0„ 

cos  ft) 
cos0„ 


p 


pdSy^dSy 


\X=const. 


(13) 


Equation  (13)  is  valid  for  o)  <  n/2.  In  the  calm  sea  limit, /7(Sx,  Sy,  0)  approaches  6(5^,  5y),  the 
double  integral  in  the  denominator  becomes  coso)/cos0n,  and  q  equals  p  regardless  of  beam 
direction. 


We  emphasize  that  we  have,  thus  far,  considered  the  interaction  of  facets  with  a  beam  point¬ 
ing  in  an  arbitrary  direction.  From  this  point  on,  however,  we  will  often  be  considering  the  situa¬ 
tion  when  the  beam  points  toward  the  receiver  (u  =  iv). 

A  plot  of  ^  as  a  function  of  slope  is  shown  in  figure  6  for  a  windspeed  of  10  m  s“^  and  a 
beam  pointing  to  zenith.  Figure  7  shows  a  plot  of  q  for  a  beam  direction  of  (80°,  270°);  for 
example,  a  receiver  due  west  of  an  ocean  scene  elevated  10°  above  the  footprint  (a  receiver  on 
the  east  coast  of  the  United  States  looking  down  at  the  Atlantic  Ocean)  with  the  wind  coming  out 
of  the  south.  Compared  to  the  distribution  of  p  whose  peak  is  at  the  origin  (figures  2,  3,  and  6), 
the  distribution  of  q  is  tilted  so  that  its  peak  no  longer  occurs  at  the  origin.  Slopes  with  positive  y 
values  are  by  and  large  responsible  for  the  interaction  in  this  case  because  they  are  tilted  toward 
the  receiver.  Those  tilted  away  from  the  receiver  (such  as  the  one  shown  in  figure  4)  are  less 
likely  to  interact  with  the  beam  {dB  will  be  smaller)  and  may  even  be  excluded  (as  the  facet  in 
figure  4  would  be  for  beam  elevations  much  lower  than  shown  there).  Figure  8  shows  the  dis¬ 
tribution  of  q  for  a  beam  direction  of  (89.75°,  270°);  for  example,  a  receiver  on  the  shoreline  of 
the  Pacific  coast  looking  due  west  with  the  wind  coming  out  of  the  north. 

As  shown  by  (13),  the  volume  in  slope  space  under  each  of  these  q  surfaces  is  unity. 
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Figure  6.  A  plot  of  q  versus  slope  for  a  windspeed  of  10  m  and  a  beam  pointing  to 
ward  the  zenith.  This  figure  is  identical  to  figure  3  apart  from  the  vertical  scale. 
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Figure  7.  A  plot  of  q  versus  slope  for  a  windspeed  of  10  m  s“^  and  a  beam  pointing  in 
the  direction  u  =  (d,  <t>)  =  (80°,  270°).  Figure  4  has  been  drawn  for  the  same  beam 
direction.  Note  that  this  distribution  is  tilted  toward  the  beam  direction  and  that  slopes 
with  large  negative  Y  values  do  not  interact  with  the  beam,  that  is,  q  is  zero  for  those 
slopes.  This  can  be  pictured  by  imagining  what  happens  in  figure  4  as  the  facet  shown 
there  is  tilted  more  and  more  toward  the  Y  direction:  the  Y  slope  becomes  more 
negative,  dB  becomes  smaller,  and  eventually  dF  contains  u  at  a  Y  slope  of  about  -0.2 
whereupon  the  facet  is  parallel  to  the  beam.  Further  increase  in  facet  tilt  exposes  the 
back  of  the  facet  to  the  beam  and  by  definition  such  facets  are  dropped  from 
consideration  and  excluded  from  interacting  with  this  particular  beam. 
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Figure  8.  A  plot  of  q  in  slope  space  for  a  windspeed  of  10  m  s~^  and  a  beam  pointing  in 
the  direction  (89.75°,  270°).  Only  facets  with  positive  y  slopes  are  allowed  to  interact 
with  this  beam. 

5.  RADIANCE  REFLECTED  BY  A  RUFFLED  FOOTPRINT 

Having  found  how  a  beam  interacts  with  a  group  of  facets,  we  will  next  show  how  those 
facets  reflect  radiance  into  a  receiver.  Figure  9  shows  the  instantaneous  position  of  a  single  facet 
whose  area  may  be  on  the  order  of  one  mm^.  The  footprint  area,  on  the  other  hand,  is  on  the 
order  of  10  m^  after  projection  normal  to  the  receiver.  Therefore,  figure  9  is  drawn  so  that  the 
facet  fails  to  fill  the  receiver  field  of  view.  The  situation  in  figure  9  is  analyzed  by  (1)  converting 
the  radiance  incident  on  a  single  facet  into  a  power,  (2)  finding  how  the  facet  reflects  that  power, 
(3)  adding  up  the  reflected  power  from  all  facets  within  the  footprint,  and  (4)  converting  the 
resultant  sum  into  a  radiance  leaving  the  footprint  for  the  receiver. 

(1)  The  source  power  arriving  at  the  facet  is  contained  in  a  beam  whose  radiance  is  Ns-  If  this 
radiance  were  received  directly  (without  being  altered  by  the  facet),  it  would  have  an  angu¬ 
lar  divergence  of  dB/R^  at  the  receiver  aperture  where  dB  is  the  facet  area  projected  normal 
to  the  receiver  and  R  is  the  distance  of  the  facet  from  the  receiver.  Let  the  area  of  the 
receiver  aperture  be  dS.  Then,  from  the  definition  of  radiance  as  the  product  of  the  beam 
power  per  unit  solid  angle  per  unit  area  normal  to  the  beam,  we  have 

dP,  =  N.^dS  =  N,^dB.  (14) 
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Figure  9.  Instantaneous  specular  reflection  by  a  single  facet.  The  radiance  of  the 
marine  sky  is  redirected  by  the  facet  into  a  receiver  looking  down  onto  the  ocean 
surface.  The  power  is  unchanged  except  for  a  reduction  due  to  the  less  than  perfect 
reflectivity  of  the  facet  at  this  particular  angle  of  incidence. 


(2)  The  facet  reflects  source  power  dPs  from  direction  Uj  to  direction  iv-  After  leaving  the  facet, 
this  power  is  multiplied  by  the  facet  reflectivity  p : 

dPr^pdP,.  (15) 

The  effect  of  the  facet  mirror  is  simply  to  alter  the  direction  of  the  power  (as  indicated  by 
the  dashed  lines  in  figure  9)  and  multiply  it  by  p . 

(3)  The  footprint,  as  a  whole,  will  contain  facets  of  all  possible  slopes.  The  power  reflected 
from  the  entire  footprint  will  be  given  by  integration  of  (15)  over  the  footprint  area  A: 

(16) 

A  A  ^ 

(4)  This  total  reflected  power  is  contained  in  solid  angle  dS/R^  and  passes  through  area 
B(Qr,  <jv>  ^  ®  Br,  the  area  of  all  facets  projected  toward  the  receiver.  Hence,  the  received 
radiance  contributing  to  a  single  pixel  is 

B-(dSIR^) 
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Substituting  (16)  foT  into  (17)  and  using  (7),  the  received  radiance  can  be  expressed  as 

^r=  llpNiqrdS^dSy 

JJn  ('«) 

u  =const 

r 

where  the  subscript  on  q  is  used  as  a  reminder  that  the  interaction  probability  density  must  be 
appli  0  a  beam  projected  toward  the  receiver.  Equation  (18)  shows  that  the  interaction  density 
q,  the  mathematical  role  of  a  weighting  factor  attached  to  the  facet  position  shown  in 

1 


<gure  10  shows  how  objects  in  the  marine  sky  contribute  to  the  radiance  of  a  pixel  in  the 
ocean.  The  conditions  are  meant  to  be  approximately  the  same  as  those  used  to  draw  figure  7, 
namely  a  receiver  elevated  by  10°  and  a  windspeed  of  10  m  s“^  Facets  of  the  various  slopes 
shown  in  figure  7  have  captured  a  broad  portion  of  the  sky  and  redirected  it  toward  the  receiver 
where  it  is  sensed  as  a  single  undivided  radiance.  The  importance  of  a  particular  slope  in  the  cap¬ 
ture  process  is  given  by  its  weighting  factor  qr  in  figure  1.^  The  full-width  half-maximum 
(FWHM)  contour  for  the  capture  is  elliptical,  rather  than  circular,  because  the  standard  devi¬ 
ations  in  the  upwind  and  crosswind  directions  are  different. 

TXvo  comments  are  offered  on  this  derivation  of  footprint  radiance.  First,  we  have  broken  the 
radiance  apart  into  a  power  and  later  reassembled  it  into  a  radiance  because  each  facet  fails  to  fill 
the  receiver  field  of  view.^  Second,  in  (17),  we  have  projected  the  facets  (rather  than  the  foot¬ 
print)  toward  the  receiver.  This  removes  the  unphysical  infinity  mentioned  in  the  introduction 
because,  although  the  footprint  projects  toward  zero  for  a  horizontal  beam,  the  facets  do  not.  The 
footprint  is  a  rough  surface,  not  a  smooth  one;  the  wind  always  tilts  some  facets  so  they  can  be 
seen  by  the  receiver  even  though  the  receiver  might  be  looking  out  in  a  direction  that  grazes 
mean  sea  level 


^  Further  study  of  the  relationship  between  figures  7  and  10  shows  that  we  have  allowed  some  slopes  to  reflect  the 
distant  sea  into  the  receiver,  but  this  is  incorrect.  Such  reflections  fall  into  the  category  of  multiple  reflections  that 
we  have  explicitly  excluded.  However,  in  the  next  section  we  will  transfer  our  attention  from  integration  over 
slopes  in  the  sea  to  integration  over  the  sky  itself  and  then  the  offending  slopes  are  conveniently  excluded  by  in¬ 
tegration  down  to,  but  not  beyond,  the  horizon. 

^  The  concept  of  radiance  invariance  requires  that  the  field  of  view  of  the  detector  be  completely  filled  by  the  distant 
surface  through  which  the  conserved  quantity  (radiance)  passes  (Preisendorfer,  1976,  vol  II,  p.  46  ff.). 

^  Please  refer  to  appendix  A  for  further  discussion. 
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Figure  10.  The  pixel  footprint  seen  by  the  receiver  of  figure  9.  Conditions  are  the  same 
as  for  figure  7.  Due  to  the  variety  of  slopes  possessed  by  all  the  facets  inside  the 
footprint,  a  broad  elliptical  portion  of  the  marine  sky  of  about  10°  in  each  direction  is 
reflected  toward  the  receiver  where  it  is  sensed  as  a  single  undivided  radiance.  The 
dotted  line  labeled  “FWHM”  indicates  those  positions  in  the  sky  that  are  reflected  into 
the  receiver  by  slopes  in  figure  7  for  which  q  has  the  value  3V2,  half  its  maximum  value 
of  7. 

6.  TRANSFORMATION  FROM  OCEAN  TO  SKY:  THE  TOLERANCE 

ELLIPSE 

Although  equation  (18)  is  the  easiest  way  to  visualize  the  geometry  of  facet  reflection 
throughout  the  footprint,  it  is  sometimes  easier  to  integrate  over  objects  in  the  sky  (rather  than 
slopes  in  the  ocean)  for  the  purpose  of  calculating  the  received  radiance.  For  a  fixed  receiver 
position,  the  location  of  the  source  radiance  in  figure  9  is  determined  by  the  slope  of  each  facet. 
Equations  can  be  derived  that  give  the  sky  location  (65,  ())j)  reflected  into  a  receiver  located  at 
(6^,  <jv)  by  a  facet  whose  slope  is  (.^,  .^).  These  equations  can  be  regarded  as  a  transformation 
from  slope  coordinates  to  sky  coordinates  with  the  receiver  position  a  parameter  of  the  trans¬ 
formation.  Hence,  in  place  of  equation  (18),  we  can  equally  well  write: 

^r=  .  (19) 

=const 
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where  we  have  introduced  the  Jacobian  J  of  the  transformation  between  facet  slope  and  sky  posi¬ 
tion.  In  appendix  A,  it  is  shown  that 


^  )  sec  (o  sin  9^  sec^  9^ 


The  geometric  meaning  of  this  Jacobian  was  originally  pointed  out  by  Cox  and  Munk  (1954) 
reference.  The  Jacobian  is  related  to  the  area  S  of  the  tolerance  ellipse:^ 

// rfs/s  =/j./rfe/0^.  ,21) 

ellipse  ^  sun 

The  word  “ellipse”  under  the  first  double  integral  means  that  the  slope  coordinates  are 
restricted  to  the  tolerance  ellipse.  The  word  “sun”  under  the  second  double  integral  means  that 
the  source  coordinates  are  restricted  to  the  solar  disk.  During  integration  over  the  solar  disk,  it  is 
even  more  convenient®  to  use  a  new  coordinate  system  X )  where  |  is  the  angular  distance 
from  the  center  of  the  disk,  and  x  is  the  azimuth  from  the  upper  limb  of  the  disk.^  Then  there 
will  be  another  Jacobian  J'  for  the  transformation  from  sky  to  disk.  It  is  shown  in  appendix  A 
that 

r._  _  sin^  (22) 

di^,x)  sin0^ 

so  that 

S  =  jjy J'd^dx  =  JJasin^  d^dx  ,  (23) 

sun  sun 

where 

secajsec^0„ 


where 


Since  the  solar  disk  is  small,  o  will  remain  almost  constant  during  the  integration  over 
and  dx  and 


S-^Gojjsm^d^dx^^GoTce^  , 


where  e  is  the  angular  radius  of  the  solar  disk  and  the  subscript  on  o  means  that  it  is  to  be 
evaluated  for  a  ray  specularly  reflected  from  the  center  of  the  sun  into  the  receiver. 

Hence,  the  product  of  J/sin0j  with  the  area  of  the  sun  is  approximately  the  area  of  the  toler¬ 
ance  ellipse. 


The  tolerance  ellipse  contains  all  those  slopes  capable  of  reflecting  a  ray  from  any  part  of  the  solar  disk  into  the 
receiver.  The  edge  of  the  tolerance  ellipse  is  given  by  equations  (A12)  and  (A13)  as  the  source  ray  negotiates 
the  periphery  of  the  solar  disk. 

8  For  the  remainder  of  this  report,  we  will  arbitrarily  switch  back  and  forth  between  integration  over  slopes,  integra¬ 
tion  over  the  scene,  and  integration  over  the  solar  disk. 

9  This  coordinate  system  is  shown  in  figure  A-5. 
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7.  SEPARATION  OF  RADIANCE  INTO  SKY  REFLECTIONS  AND  SUN 

GLINTS 

We  now  return  to  the  calculation  of  the  radiance  received  from  the  fluctuating  surface  of  the 
wind-ruffled  open  ocean.  Inserting  equation  (11)  for  9;.  and  equation  (20)  for  7  into  equation 
(19),  we  arrive  at  the  following  general  expression  for  the  radiance  received  from  the  footprint; 

^  u ,  -const 

Here,  Br  is  given  by  (12)  with  the  beam  pointing  toward  the  receiver  (that  is,  (12)  with 
u  =  Ur)  and  p{Sx,  Sy,  W)  is  given  by  (2).  In  equation  (26),  the  beam  from  the  source  is  specularly 
reflected  to  the  receiver;  therefore,  ca  is  given  by  equation  (A15),  0„  is  given  by  equation  (A16) 
or  (A17),  Sx  is  given  by  equation  (A12),  and  Sy  is  given  by  equation  (A13).  Nr  is  independent  of 
source  coordinates  and  facet  slope  but  still  depends  on  sky  radiance,  sea  water  reflectivity,  wind 
velocity,  and  receiver  position.  Equation  (26)  is  the  integral  equation  for  the  receiver  radiance 
we  have  been  seeking.  It  is  an  exact  mathematical  consequence  of  the  Cox-Munk-Plass  model 
for  the  wind-ruffled  sea  surface. 

Next,  we  separate  the  received  radiance  Nr  into  two  parts:  N;--*^due  to  the  sky  with  the  sun 
removed,  and  due  to  the  sun  alone: 


N,  =  Nfy  + 

(27) 

(28) 

sfy 

u -const 

JV;""=  ]\pN,qrJde,d^, 

(29) 

sun 

u^=const 

Equation  (28)  gives  the  radiance  received  from  the  entire  sky  dome  (with  the  sun  removed) 
after  reflection  by  the  wind-roughened  surface  of  the  sea.  This  expression  cannot  be  simplified 
further  without  adopting  a  model  for  <})j),  the  sky  radiance  before  reflection.  However,  if 
the  sky  radiance  prior  to  reflection  is  independent  of  azimuth,  the  reflected  sky  radiance  can  be 
written  as 

1  j4  231 

^  "r  0  0 

where  it  is  understood  that  the  receiver  coordinates  must  be  held  constant  during  these  integra¬ 
tions. 

Equation  (29)  gives  the  sun  glint,  the  radiance  received  from  the  sun  alone  after  reflection  by 
the  ruffled  sea  surface.  Here,  the  source  radiance  Ns  is  equal  to  the  solar  radiance  Nq.  Since  solar 
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radiance  is  constant  throughout  the  solar  disk  (the  sun  is  a  Lambertian  source),  c^n  be 
brought  outside  the  integral  and  we  have  the  following  equivalent  expressions  for  the  ratio  of 
glint  radiance  leaving  the  footprint  to  solar  radiance  arriving  at  the  footprint:*^ 


for  integration  over  the  ocean, 


j^sun 

r 


®  ellipse 

u ,  -const 


i_d 

4B, 


|j  p  sec'* 

sun 

u,  =const 


d„  p  sin  Os  de,  dip, 


(31) 


(32) 


for  integration  over  the  sky,  and 


^sun 


for  integration  over  the  sun. 


iJP9r  asxn^d^dx 

disk 

Uf=const 


(33) 


In  the  calm  sea  limit,  q  is  equal  to  p,p  approaches  a  delta  function  selecting  zero  slope,  and 
we  see  from  (31)  that 


as  expected. 


(34) 


Momentarily  setting  the  reflectivity  of  the  sea  surface  to  100%  for  the  purposes  of  illustra¬ 
tion,  we  see  that  the  radiance  ratio  in  (31)  is  given  by  the  volume  of  the  interaction  probability 
density  over  the  tolerance  ellipse.  This  is  illustrated  in  figure  li.  Figure  11  is  an  expanded  view 
of  figure  7,  which  was  drawn  for  the  receiver  looking  straight  toward  the  sun  (more  precisely, 
for  an  incident  azimuth  of  90°  and  a  reflected  azimuth  of  270°).  There  is  one  elliptical  column 
in  figure  11  for  solar  zenith  angles  of  70°,  80°,  and  90°,  and  the  location  of  each  column  within 
the  figure  is  determined  by  the  solar  position.  Recalling  that  the  receiver  beam  zenith  angle  is 
80°  for  this  figure,  a  solar  zenith  angle  that  is  also  80°  requires  facets  with  slopes  close  to  zero 
for  specular  reflection.  Hence,  the  column  for  the  reflection  at  80°  is  centered  at  the  origin  of 
slope  space.  Solar  zenith  angles  smaller  than  80°  (sun  higher)  require  facets  with  positive  y  slope 
for  specular  reflection;  solar  zenith  angles  larger  than  80°  (sun  lower)  require  facets  with  nega¬ 
tive  y  slope  for  specular  reflection.  The  ratio  of  glint  radiance  leaving  the  footprint  to  solar 
radiance  arriving  at  the  footprint  is  given^^  by  the  volume  of  each  column. 


Recall  that  we  are  neglecting  the  transmission  of  the  atmosphere  so  the  radiance  at  the  top  of  the  atmosphere  and 
at  the  footpint  are  the  same. 

Ignoring  the  effects  of  reflectivity. 
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y  SLOPE  -0.05 
(CROSSWIND) 


-0.05  ^  SLOPE 

-0.1" -0.1  (UPWIND) 


Figure  11.  An  expanded  view  of  figure  7  with  glint  columns  included.  The  base  of  each 
column  is  the  tolerance  ellipse  and  the  volume  (neglecting  reflectivity  losses)  of  each 
column  is  the  ratio  of  glint  radiance  leaving  the  footprint  to  solar  radiance  arriving  at  the 
footprint.  The  receiver  is  fixed  at  (80°,  270°).  The  solar  position  for  each  column  is  (dg, 
9Cr)  where  6s  is  given  at  the  base  of  each  column  and  the  solar  azimuth  has  been  cho¬ 
sen  so  that  the  receiver  is  looking  directly  along  the  center  of  the  glint  pattern. 

Figure  12  shows  the  same  q  surface  as  that  shown  in  figure  8.  In  figures  8  and  12,  the  zenith 
angle  of  the  beam  leaving  the  footprint  for  the  receiver  is  89.75°  so  the  received  beam  almost 
grazes  the  ocean  surface.  In  figure  12,  we  have  added  crescent  caps,  which  are  the  tops  of  each 
glint  column  for  solar  zenith  angles  stepping  in  2°  increments  from  78°  to  90°.  As  the  sun 
approaches  the  horizon,  the  width  of  the  tolerance  ellipse  in  the  Y  direction  remains  constant,  but 
the  width  in  the  direction  approaches  infinity.  At  the  same  time  the  position  of  the  glint 
column  approaches  the  origin  of  slof>e  space  where  the  height  of  the  q  surface  becomes  very 
small  regardless  of  slope.  Hence,  even  though  the  area  of  the  ellipse  becomes  infinite  (as  would 
be  expected  from  equation  (20)  for  o)  equal  to  Ji/2)  the  glint  radiance  ratio  given  by  (31)  remains 
finite  for  a  horizontal  view  of  the  setting  sun. 


12  This  can  be  seen  by  imagining  a  receiver  on  the  shore  looking  directly  into  a  sun  whose  center  is  exactly  on  the  far 
horizon.  Only  the  upper  half  of  the  sun’s  disk  would  be  visible.  A  facet  with  zero  slope  would  be  required  to 
reflect  the  very  top  of  this  half-disk  into  the  receiver,  but  facets  with  infinite  slope  would  be  required  to  reflect  into 
the  receiver  those  sides  of  the  half-disk  that  touch  the  horizon. 
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Figure  12.  Figure  8  showing  the  caps  of  a  line  of  glint  columns.  The  zenith  angle  of  the 
solar  center  progresses  from  78°  to  90°  in  2°  steps  at  a  fixed  azimuth  of  90° .  The  last 
column  occurs  when  both  sun  and  receiver  are  very  close  to  the  horizon.  It  has  a  base 
that  is  very  large  (approaching  oo  in  the  positive  and  negative  x  directions),  but  it  still 
has  a  finite  volume  since  q  approaches  zero  in  those  directions.  The  approximation 
associated  with  equation  (35)  amounts  to  the  assumption  that  each  coiumn  has  a  flat 
top. 

The  behavior  of  equations  (31)  through  (33)  for  a  receiver  looking  straight  toward  a  setting 
sun  is  given  in  figures  13  through  15  by  the  curves  labeled  “Exact.”  Exact  glint  images  are  also 
shown  in  figures  16  and  17  for  a  receiver  pointed  into  a  sun  elevated  at  2°  over  moderate  and 
calm  seas,  respectively.  Note  that  the  vertical  scale  in  these  figures  never  exceeds  2%,  which 
may  be  surprising  in  view  of  the  calm  sea  limit  set  forth  in  (34),  especially  for  a  sun  setting  on 
the  horizon  for  which  p  approaches  100%.  The  exact  calculations  produce  a  smaller  value  than 
100%  because  most  of  the  slopes  in  a  footprint  deflect  solar  rays  away  from  the  receiver  toward 
other  parts  of  the  sky  even  at  very  low  windspeeds.  In  geometric  terms,  the  volume  of  qr  is  al¬ 
ways  much  larger  than  the  volume  of  a  typical  glint  column.  Put  another  way,  the  real  sea  always 
contains  an  appreciable  number  of  slopes  in  excess  of  the  solar  radius. 

We  close  this  section  with  a  few  remarks  regarding  the  relative  strengths  of  the  reflected  sky 
and  sun  glint  radiance  in  the  visible  region.  The  first  originates  in  a  relatively  weak  source, 
namely  scattered  solar  radiation.  The  second  originates  in  a  very  strong  source,  namely  the  sun 
itself.  However,  the  integral  in  (28)  for  the  sky  covers  the  entire  hemisphere  of  the  sky  (2ji  ^r), 
whereas  the  integral  in  (29)  for  the  glint  covers  only  the  solar  disk  (6,77  x  10"^  sr).  The  amount 
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of  sky  actually  contributing  to  the  integral  in  (28)  will  depend  on  how  rapidly  q  falls  off  with 
slope.  For  calm  conditions,  a  relatively  small  region  of  the  sky  near  the  specularly  reflected  (zero 
slope)  receiver  location  will  contribute  to  the  integral.  For  rough  conditions,  many  more  slopes 
will  be  present  and  a  much  larger  region  of  the  sky,  more  like  that  shown  in  figure  10,  will  con¬ 
tribute  to  the  integral.  The  point  remains  that  the  glint  and  reflected  sky  contributions  can  be 
comparable  in  the  visible  part  of  the  spectrum  even  though  the  radiance  of  the  sun  is  many 
orders  of  magnitude  greater  than  the  radiance  of  the  sky.  For  example,  on  a  clear  day  when  the 
sun  shines  with  maximum  radiance,  the  irradiance  in  the  visible  part  of  the  spectrum  at  the  sur¬ 
face  of  the  earth  due  to  the  sun  is  only  about  10  times  the  irradiance  due  to  the  sky.^^  On  cloudy 
days  (or  in  the  infrared  part  of  the  spectrum  where  solar  radiance  drops  by  an  order  of  magni¬ 
tude),  the  ratio  is  even  smaller. 


Figure  13.  Glint  ratio  versus  solar  zenith  angle  for  a  horizontal  view  of  a  sunset  over  a 
moderate  sea.  Parameters  are  the  same  as  for  figure  12  with  p  set  to  100%.  The  solid 
line  labeled  “Exact”  gives  the  volume  of  the  glint  columns  whose  caps  are  shown  in  fig¬ 
ure  12.  The  curves  labeled  “Flattop”  and  “Cox-Munk”  show  the  behavior  of  equations 
(35)  and  Cox  and  Munk  (1954  ((9)),  respectively. 


This  fact  was  pointed  out  to  us  by  Isador  Halberstam,  Hughes  STX,  Lexington,  MA.  Thank  you,  Izzy. 
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Figure  14.  Ratio  of  sun  glint  radiance  leaving  the  footprint  to  solar  radiance  arriving  at 
the  footprint  as  a  function  of  the  zenith  angle  of  the  beam  leaving  the  footprint  for  the 
receiver.  Conditions:  (65.  <Ps)  =  (88°,  90°),  <l>r  =  270P,  W=  10  m  s"^  and  p  =  100%. 


Figure  15.  Ratio  of  glint  radiance  leaving  the  footprint  to  solar  radiance  arriving  at  the 
footprint  as  a  function  of  windspeed.  In  this  figure  {Og,  <i>s)  =  (8&‘,  90°),  (Or,  <pr)  =  (88°, 
270°),  and  p  =  100%. 
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Figure  16.  Glint  pattern  for  a  sun  elevated  by  2°  over  a  moderate  sea.  The  sun  is  to  the 
right  of  the  diagram,  near  the  positive  Y  axis,  with  an  azimuth  of  90°.  The  receiver  is  to 
the  left  of  the  diagram,  near  the  negative  Y  axis,  looking  directly  toward  the  sun.  The 
windspeed  is  10  m  s~^.  The  reflectivity  is  100%. 

z 


ZENITH  ANGLE  (deg) 
OF  RECEIVER  BEAM 


AZIMUTH  (deg)  OF 
RECEIVER  BEAM 


Figure  1 7.  Glint  pattern  for  a  sun  elevated  by  2°  over  a  calm  sea.  Except  for  a  wind- 
speed  of  1  m  conditions  are  the  same  as  for  figure  16.  Note  the  change  in  vertical 
scale  compared  to  the  previous  figure. 
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8.  THE  “FLATTOP”  APPROXIMATION  FOR  THE  GLINT  RADIANCE 


Figures  13  through  15  also  contain  curves  marked  “Flattop”.  These  are  an  approximation, 
which  will  now  be  derived,  for  the  exact  glint  expressions.  The  idea  behind  the  “Flattop” 
approximation  was  used  earlier  without  comment  to  obtain  (25)  as  an  approximation  for  the 
tolerance  ellipse;  the  smalt  size  of  the  solar  disk  allows  an  integrand,  this  time  the  one  in  (31),  to 
be  brought  outside  an  integral.  Then,  the  ratio  of  glint  radiance  leaving  the  footprint  to  solar 
radiance  arriving  is  given  approximately  by 


j^sun 

-TJ—  “  { Wr }o  IJ y={pqr}o  -^ 

°  ellipse 

n,=comt 


for  integration  over  the  ocean,  or 
j,sun 

No  ^45, 


1  A  4 


o  sun 
u=const 


sec^  6„p 

[4  5, 


1  ^ 

sec  0„  p 


\\sm^d^dx 


[4  5, 


for  integration  over  the  sky. 


o  disk 
u,,  ^const 

o 


(35) 


(36) 


The  meaning  of  (35)  can  be  appreciated  by  looking  at  figures  11  and  12.  This  approximation 
(ignoring  the  effects  of  reflectivity)  amounts  to  replacement  of  the  sloping  column  cap  by  a  flat- 
top  whose  height  is  the  value  of  q  at  the  column  center. 

This  “Flattop”  approximation  fails  when  the  integrand  varies  appreciably  during  integration 
over  the  solar  disk,  a  situation  which  occurs  when  both  the  sun  and  the  receiver  are  near  the  hori¬ 
zon.  Then,  the  problem  with  the  approximation  can  easily  be  seen  in  figure  12.  As  the  receiver 
approaches  the  horizon  with  the  sun  also  low  in  the  sky,  the  width  of  the  tolerance  ellipse 
approaches  infinity  in  the  A"  direction  giving  rise  to  a  substantial  overestimate  of  the  column 
volume  when  using  the  Flattop;  the  actual  column  height  used  in  the  exact  calculation  varies 
appreciably  across  the  elliptical  base  approaching  zero  in  the  wings.  As  a  result,  when  both 
source  and  receiver  are  near  the  horizon,  the  real  volume,  and  the  radiance  it  represents,  is  less 
than  the  approximate  value  given  by  (35).  These  trends  are  evident  in  figures  13  and  14- 

It  is  obvious  from  a  comparison  of  figures  2  and  3,  that  the  “Flattop”  approximation 
improves  as  the  windspeed  increases.  Although  we  have  not  made  a  study  of  the  detailed 

The  column  center  is  the  specular  reflection  slope  for  the  center  of  the  sun. 
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conditions  under  which  this  approximation  is  valid,  we  will  note  that  the  approximation  remains 
finite  even  when  evaluated  for  infinite^^  slopes  where  6„  approaches  Ji/2.  This  can  been  under¬ 
stood  by  referring  to  (36)  and  noting  that:  {\)AIBr  remains  finite  for  receivers  on  the  horizon,^^ 
(2)  sec^Ofl  approaches  infinity  like  where  R  is  the  radius  in  slope  space  to  the  center  of  the 
column,  and  (3) p  falls  off  for  large  slopes  like  cxp{-R^lcs^)  where  o  is  the  rms  average  of 
upwind  and  crosswind  deviations. 


9.  THE  “COX-MUNK”  APPROXIMATION  FOR  THE  GLINT  RADIANCE 


The  “Cox-Munk”  approximation  can  be  introduced  by  noting  that  facets  tilted  by  less  than 
the  beam  elevation  will  always  present  their  front  surface  to  the  beam.  If  the  wind  velocity  is 
low  enough  that  area  A  is  almost  entirely  occupied  by  such  facets,  the  restriction  on  the  angle  of 
incidence  in  (12)  and  (13)  may  be  dropped  without  appreciable  error,  and  we  may  write  (see 
appendix  A)  that 


5 

A 


for  an  arbitrary  beam,  and 


// 

u -const 


COSO) 

cos6„ 


pdS^dSy  =  COS0 


—  =  cos0^ 
A  ‘ 


(37) 


(38) 


for  a  beam  pointing  toward  the  receiver. 


Substituting  (38)  into  (36),  we  obtain 


N 


sun 


1 


Ng-Tte^  14  J 

and  since  No*Tce^  is  equal  to  the  solar  irradianc«  Hq,  we  also  have 

“  secOr  p  sec"^  d„ 


N 


sun 


(39) 


(40) 


for  the  ratio  of  glint  radiance  leaving  the  ocean  surface  to  the  solar  irradiance  arriving  at  the 
ocean  surface.  Equation  (40)  is  identical  to  the  equation  used  by  Cox  and  Munk  (1954  (9)). 


The  “Cox-Munk”  approximation  is  valid  under  the  same  conditions  for  which  (38)  is  valid. 
A  detailed  study  of  those  conditions  is  given  in  appendix  A.  In  contrast  to  the  “Flattop”  approxi¬ 
mation,  the  “Cox-Munk”  approximation  improves  as  the  windspeed  decreases;  for  a  receiver 
beam  elevated  by  10°  or  more,  the  error  is  less  than  10%  for  windspeeds  below  11  m  s~^  and 
less  than  1%  for  windspeeds  below  SVi  m  s~^ 


15 


Cf.  figures  A-1  through  A-4  and  the  discussion  are  in  appendix  A. 
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10.  CONCLUSION 

Saunders  (1968),  in  the  spirit  of  the  first  class  of  solutions  mentioned  in  the  introduction, 
pointed  out  that  multiple  retlections  and  shadowing  should  be  included  in  any  theory  of  solar 
glints.  Multiple  reflections  occur  when  the  source  ray  bounces  off  of  several  facets  before  reach¬ 
ing  the  receiver.  Shadowing  refers  to  the  fact  that  slopes  on  the  back  sides  of  the  waves  and  deep 
in  the  troughs  between  waves  are  hidden  from  view.  Each  of  these  effects  becomes  more  impor¬ 
tant  when  the  source  or  receiver  approaches  the  horizon.  (Cox  and  Munk  did  not  include  either 
effect  but  were  careful  to  take  measurements  when  both  source  and  receiver  were  near  the  zenith 
to  minimize  the  contributions  of  each  effect.)  Saunders  established  upper  and  lower  limits  to  the 
observed  radiance  based  on  multiple  scattering  and  he  included  shadowing  by  introducing  a 
shadowing  factor.  His  shadowing  factor,  which  multiplies  the  equation  (Cox  and  Munk,  1954 
(9)),  clamps  the  Cox-Munk  radiance  to  a  finite  value  for  horizontal  observation  of  solar  glints. 
However,  we  do  not  find  that  the  use  of  Saunders’  shadowing  factor  produces  the  proper 
radiance  in  the  calm  sea  limit. 

In  contrast  to  the  approach  of  Saunders  (1968)  to  removing  the  unphysical  infinity,  we  have 
followed  a  suggestion  by  Plass  et  al.  (1975).  Neglecting  multiple  scattering  and  shadowing  just 
as  Cox  and  Munk  did,  we  have  obtained  a  new  integral  version  of  Cox  and  Munk  (1954  (9)), 
which  is  correct  for  all  geometries.  The  integral  version,  equation  (32),  produces  a  finite  glint 
radiance  reflected  horizontally  by  the  ocean  surface  even  when  the  sun  is  also  on  the  horizon. 
Furthermore,  it  gives  the  proper  radiance  in  the  calm  sea  limit. 

When  the  reflected  beam  is  elevated  by  more  than  about  10®  over  a  moderate  sea,  an  approx¬ 
imation  may  be  employed  that  reduces  the  exact  formulation  we  have  derived  to  (40),  the  form 
used  by  Cox  and  Munk.  We,  therefore,  see  that  the  resolution  of  the  unphysical  infinity 
motivating  this  report  lies  in  the  third  solution  class  of  the  introduction.  The  Cox-Munk  radiance 
equation,  equation  (40),  is  an  approximation  valid  for  radiance  away  from  the  horizon.  Near  the 
horizon,  integral  equation  (32)  should  be  used  instead  of  approximate  equation  (40). 

General  equation  (26),  valid  for  all  viewing  geometries,  will  provide  the  proper  foundation 
for  an  improved  model  of  ocean  radiance  incorporating  multiple  scattering  and  shadowing. 
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12.  SYMBOLS 

Notes:  Column  I:  Cox-Munk  (1954).  Column  II:  this  report,  text.  Column  III:  this  report, 
FORTRAN  source  code.  All  symbols  connected  with  radiation  (for  example  radiance,  irradiance, 
and  reflectivity)  stand  for  spectral  quantities  valid  at  a  single  (unspecified)  wave  number. 


I 

II 

III 

Latin 

A 

A 

Area  of  pixel  footprint. 

a 

a 

X  component  of  unit  vector.  Upwind. 

B 

B 

Area  of  ergodic  cap  projected  toward  beam. 

b 

b 

Y  component  of  unit  vector.  Crosswind. 

c 

c 

Z  component  of  unit  vector.  Zenith. 

dA 

Area  of  single  facet  projected  onto  horizon. 

dB 

Area  of  single  facet  projected  toward  beam. 

H 

Irradiance  (W m~^  [cm~^]“^). 

H 

Ho 

Solar  irradiance  arriving  at  footprint. 

J 

Jacobian  of  transformation,  slopes  to  sky. 

r 

Jacobian  of  transformation,  sky  to  disk. 

N 

N 

Radiance  (lY  sr“^  [cm“^]“^). 

Nr 

Reflected  radiance. 

Nr^'^y 

Reflected  radiance  from  sky. 

N 

j^^sun 

Reflected  radiance  from  sun. 

Emitted  radiance  from  sea. 

No 

No 

Incident  radiance  of  sun. 

Ns 

Incident  radiance  of  source  (sky  or  sun). 

o 

o 

Subscript  referring  to  the  center  of  the  sun. 

p 

Probability  of  facet  occurrence. 

p 

Beam  power  (W). 

P 

p 

P 

Probability  density  of  facet  occurrence. 

Q 

Probability  of  facet-beam  interaction. 

q 

q 

Probability  density  of  facet-beam  interaction. 

R 

Radius  in  slope  space. 

R 

Distance  from  facet  to  receiver. 

At 

S 

s 

Area  of  tolerance  ellipse. 

Zx 

Sx 

Sx 

Facet  slope  in  x  direction. 

Sy 

Sy 

Facet  slope  in  y  direction. 

U 

Unit  vector. 

Unit  vector  normal  to  facet. 

Ur 

Unit  vector  pointing  to  receiver. 

Us 

Unit  vector  pointing  to  source. 

W 

w 

Wind  speed. 

X 

Coordinate  axis  pointing  upwind. 
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Y 

Coordinate  axis  pointing  crosswind. 

Z 

Coordinate  axis  pointing  to  zenith. 

Greek 

a 

Azimuth  of  steepest  ascent. 

E 

E 

Epsilon 

Angular  radius  of  solar  disk. 

V 

V 

Wave  number  (cm“^). 

Solar  elevation. 

<t> 

P 

Azimuth  of  u. 

4^n 

Pn 

Azimuth  of  u„. 

<tv 

Pr 

Azimuth  of  iv- 

<j)j 

Ps 

Azimuth  of  Uj. 

e 

T 

2^nith  angle  of  u. 

0n 

Tn 

Zenith  angle  of  u„.  Tilt. 

fA 

Br 

Tr 

Zenith  angle  of  Up 

es 

Ts 

Zenith  angle  of  u^. 

P 

P 

Rho 

Reflectivity  of  sea  water. 

O 

Sigma 

Area  ratio:  (tolerance  ellipse)/(solar  disk), 

Oc 

Oc 

St 

Standard  deviation  of  Sy. 

Ou 

Ou 

Su 

Standard  deviation  of 

U) 

(0 

Omega 

Angle  of  incidence,  angle  of  reflection. 

CO 

Angle  between  beam  and  facet  normal. 
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APPENDIX  A 

MATHEMATICAL  DETAILS  OF  THE  COX-MUNK-PLASS  MODEL 
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CARTESIAN  COORDINATES 


The  coordinate  system  for  this  report  is  shown  in  figure  1.^  We  have  defined  unit  vectors 

u  =  (e,(p)  =  (a,b,c)  (Al) 

Here,  0  is  the  zenith  angle  of  the  vector,  and  <}>  is  the  azimuth  of  the  vector.  The  X,  Y,  and  Z 
components  of  the  vector  are  a,  b,  and  c  respectively.  Each  of  the  vectors  points  away  from  the 
origin.  Azimuths  are  considered  positive  when  measured  counterclockwise  looking  along  the 
nadir. 


By  resolving  polar  components  into  vertical  and  horizontal  components,  and  then  by  resolv¬ 
ing  the  horizontal  component  into  its  X  and  Y  components,  we  obtain 

a  =  sm6cos(f> 

b  =  sm  Osin  0  (A2) 

c  =  cosO 


for  each  vector.  Subscripts  attached  to  a  unit  vector  are  used  to  refer  to  the  three  vectors  involved 
in  specular  reflection;  the  subscript  will  equal  s  when  referring  to  the  source,  n  when  referring  to 
the  facet  normal,  and  r  when  referring  to  the  receiver. 

It  will  be  convenient  to  have  expressions  for  the  facet  orientation  in  terms  of  =5  z/^x,  the 

facet  slope  in  the  x  direction,  and  Sy  =  ^  z/^y,  the  facet  slope  in  the  y  direction.  These  expres¬ 
sions  can  be  obtained  by  noting  that  any  vector  r  =  (x,  y,  2)  within  the  facet  is  perpendicular  to 
the  facet  normal: 

vu„=xa„+yb„-i-zc„^0  ,  (A3) 

where 

al+b^  +cl=\.  (A4) 

n  n  n 


Equation  (A3)  sets  the  direction  of  the  normal,  and  equation  (A4)  sets  its  length.  By  evaluat¬ 
ing  equation  (A3)  in  the  X-Z  and  Y-Z  planes,  the  slopes  are  found  to  be 


(A5) 


Inserting  (A5)  into  (A4)  and  using  the  Z  component  of  (A2),  we  also  have 


+  sj  =tan^0„. 


(A6) 


*  Figures  1  through  17  and  equations  (1)  through  (40)  appear  in  the  main  body  of  this  report. 
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PROJECTED  AREA  OF  THE  ERGODIC  CAP 


If  all  the  facets,  like  the  one  5>hown  in  figure  4,  are  smoothly  joined  together,  they  form  a 
surface  that  Preiscndorfer  (1976)  has  called  the  ergodic  cap.  B  represents  the  area  of  the  ergodic 
cap  projected  in  the  arbitrary  direction  u.  That  portion  of  B  produced  by  a  single  facet  is  shown 
in  figure  4,  and  an  expression  for  BIA  is  given  by  equation  (12). 


We  will  now  present  an  approximate  expression  for  5/A,  which  is  valid  at  low  windspeeds 
for  beams  near  the  zenith.  For  low  windspeeds,  the  occurrence  probability  density  p  is  narrowly 
peaked  about  zero  facet  tilt;  p  will  approach  a  two-dimensional  delta  function.  Furthermore,  very 
large  slopes  are  required  for  facets  that  present  their  backs  to  beams  near  the  zenith,  and  such 
large  slopes  are  very  infrequent.  Hence,  for  calm  enough  conditions,  very  few  facets  will  occur 
with  their  backs  to  any  given  beam,  especially  those  beams  near  the  zenith,  and  the  restriction  on 
the  integral  in  (12)  that  a)sjt/2  can  be  dropped  without  appreciable  error.  In  this  case,  we  have 


B 

A 


IJ 

u=const 


COSO) 

COS0„ 


pdS^dSy. 


The  trigonometric  factors  in  the  integrand  of  (A7)  are 


(A7) 


COSfl) 

COS0„ 


u 


n 


COS  Or 


_  aa„  +  bb„  +  cc„ 


—  —qS^  ~  bSy  +  c . 


(A8) 


When  (A8)  is  inserted  into  (A7),  the  beam  factors  (those  involving  a,  b,  and  c)  can  be 
brought  outside  the  integral  since  they  do  not  depend  on  the  facet  slope.  Using  the  facts  that 
(1)  p  is  even,  (2)  Sx  and  Sy  are  both  odd,  and  (3)  p  is  normalized  to  unity,  the  first  two  terms  in 
(A8)  vanish  altogether  and  the  last  term  in  (A8)  survives  integration  intact.  Hence, 

COSO)  . 

- pdSxdSy=c  =  cos6  (p 

cos0„  ^  ^ 


u-const 


for  the  approximate  area  of  the  ergodic  cap  projected  in  the  arbitrary  direction  u. 

Figure  A-1  shows  how  A/B,  the  reciprocal  projected  area  of  the  ergodic  cap,  behaves  when 
projected  anywhere  in  the  sky  near  the  horizon.  In  this  region,  the  approximation  just  derived 
begins  to  fail.  For  zenith  angles  greater  than  about  80°,  A/B  departs  from  sec0  by  an  amount  that 
depends  on  windspeed.  The  higher  the  windsjseed,  the  larger  the  deviation  from  sec0.  Physically, 
this  is  because  the  wind  tilts  the  facets  away  from  their  calm  horizontal  orientation  in  all  direc¬ 
tions.  Some  are  tilted  away  from  a  beam  that  grazes  the  mean  surface  of  the  sea;  these  facets 
cannot  reflect  and  are  ignored  in  the  integral.  Some  facets  are  tilted  toward  the  grazing  beam; 
these  facets  project  more  area  toward  the  beam,  are  more  effective  reflectors  than  if  they  were 
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horizontal,  and  contribute  more  to  the  integral.  The  net  result  is  that  a  wind-ruffled  ocean  surface 
projects  a  finite  area  (rather  than  a  zero  area)  horizontally.  This  fact  is  responsible  for  many 
simple  observations  such  as  the  very  existence  of  the  visual  horizon  (Cox  and  Munk,  1955)  and 
the  absence  of  reflections  of  low-lying  distant  land  masses  (Minnaert,  1954). 


WIND 


Figure  A-1.  The  reciprocal  area  of  the  ergodic  cap  seen  by  a  beam  as  a  function  of  that 
beam’s  zenith  angle.  For  zenith  angles  less  than  about  80°  the  secant  of  the  zenith 
angle  may  be  used  for  this  reciprocal  area  to  a  fair  approximation.  For  beams  closer  to 
the  horizon,  however,  the  true  values  deviate  from  the  seed  approximation,  as  shown  in 
this  figure,  with  larger  deviations  at  larger  windspeeds.  The  beam  azimuth  is  0°. 

Figure  A-2  gives  the  same  information  as  figure  A-1  in  a  different  form:  it  shows  how  the 
ratio  ofB  toA»cos6  behaves  near  the  horizon  and  is  useful  in  estimating  the  accuracy  of  (A9) 
for  various  zenith  angles  and  windspeeds. 

Figures  A-1  and  A-2  have  been  for  a  beam  azimuth  of  0°.  Figures  A-3  and  A-4  show  the 
same  results  for  a  beam  azimuth  of  90°.  For  a  given  windspeed,  results  for  all  other  azimuths  are 
intermediate  to  the  results  given  by  these  two  cases. 

Finally,  we  note  that  for  beams  with  zenith  angles  less  than  80°  and  any  azimuth  whatsoever, 
the  error  involved  in  the  approximation  given  by  (A9)  is  less  than  1%  for  windspeeds  below  SVz 
m  s~^  and  less  than  10%  for  windspeeds  below  11  m  s~^ 
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70  75  80  85  90 

ZENITH  ANGLE  OF  BEAM  6  (deg) 


Figure  A-2.  The  ratio  between  the  area  of  the  ergodic  cap  and  A»cosd  as  a  function  of 
the  beam  zenith  angle  S.  The  more  calm  the  conditions  the  closer  to  the  horizon  the 
A*cosB  approximation  can  be  used.  The  beam  azimuth  is  0°. 


WIND 


Figure  A-3.  The  same  as  figure  A-1  except  for  a  beam  azimuth  of  90° .  All  other  azi¬ 
muths  display  deviations  intermediate  between  these  two  cases  of  0°  and  90°  azimuth. 
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f 


70  75  80  85  90 

ZENITH  ANGLE  OF  BEAM  0  (deg) 


Figure  A-4.  The  same  as  figure  A-2  except  for  a  beam  azimuth  of  90° . 


APPROXIMATE  EXPRESSION  FOR  THE  INTERACTION 
PROBABILITY  DENSITY 


When  the  beam  vector  u  points  toward  a  receiver  well  above  the  horizon  and  the  wind 
velocity  is  small,  B/A^cosQr-  Then  the  interaction  probability  density  in  equation  (11)  can  be 
approximated  by 


COSW 

- P' 

cos  6;.  cos 


(AlO) 


FACET  GEOMETRY  FOR  SPECULAR  REFLECTION 

For  fixed  positions  of  the  source  and  receiver,  the  law  of  reflection 

u„  =  (uj  +  u^)/(2cosa})  (All) 


will  determine  the  facet  position  for  a  specular  reflection  between  source  and  receiver. 
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We  require  equations  for  the  geometry  of  the  facet  (slopes,  angle  of  incidence,  and  tilt)  in 
terms  of  the  source  and  receiver  coordinates.  To  find  Sx  and  Sy,  we  can  write 


(sin  cos^5  +  sin  6^  cos^^ ) 
(cos6^  +  COS0;. ) 


(A12) 


(bs+i>r) 

Cn 

(sin  0^  sin  +  sin  0^.  sin  0^ ) 
{cosOg  +  cosdr ) 


(A13) 


The  angle  of  incidence  can  be  found  by  squaring  equation  (All)  to  obtain 

2C0S^  Ct)  =  l  +  Uj  *0;..  (A14) 

Substituting  Cartesian  coordinates  for  source  and  receiver  gives  the  following  expression  for  the 
angle  of  incidence  in  terms  of  source  and  receiver  positions: 

2cos^  (O  =  l  +  sin0^  sin0,.cos(0^  -0^)  +  cos0^cos0;.  .  (A15) 

Equations  (A12)  and  (A13)  can  be  substituted  directly  into  equation  (A6)  to  derive  this 
expression  for  the  tilt 

2  _  sin^  0j  +  sin^  6^  +  2sin0j sin0^  cos(0j  -  0r ) 

tan  0„ - - - — -  ,  (A16) 

(COS0^  +  COS0;. ) 

and  a  second  expression  for  the  tilt  can  be  found  directly  from  the  z  component  of  equation 
(All); 

COS0„  =  (cos0j  +  C0S6y)/2C0S(0  .  (A17) 

We  reemphasize  that  equations  (All)  through  (A17)  apply  only  in  case  a  ray  from  the  source 
is  reflected  by  the  facet  toward  the  receiver. 
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THE  JACOBIAN  FOR  THE  TRANSFORMATION  BETWEEN  OCEAN 

AND  SKY 


Equations  (A12)  and  (A13)  formally  express  a  mathematical  transformation  between  source 
coordinates  in  the  sky  and  slope  coordinates  in  the  sea,  if  the  receiver  position  is  regarded  as  a 
parameter  of  the  transformation.  The  Jacobian  of  the  transformation  is 


dSx/dds  dSx/d((>s 

Calculation  of  the  lal  derivatives  within  the  determinant  gives 

^2  -  cosdj  cosips  +  sin  OgSjc  +  sin  9^  sin  (p^ 

-  cos  dg  sin  (p^  +  sin  0^5^  -  sin  9^  cos  (p^ 

where  h  =  cosds  +  cos0r.  When  the  determinant  i-.  evaluated  we  find 


(A18) 


(A19) 


sin0. 


=  COS0J  -  sin  9g  {sin  0^Sy  +  cos^^Sj^ }  ■ 


(A20) 


Substituting  equations  (A12)  and  (A13)  for  the  slopes,  the  expression  between  braces 
becomes 


(-(sin  9r  cos(0s  ~  0;. )  +  sin  0^ )//?} . 


(A21) 


When  (A21)  is  multiplied  by  -sin0j  and  recast  by  means  of  (A16),  we  have  for  the  second 
term  on  the  right-hand  side  of  (A20) 

-sin 0j {  }  =  (h^  tan^  0„  +  sin^  0^  -  sin^  9^ )  jlh .  (A22) 

Using  h  -  cos0r  for  the  first  term  on  the  right-hand  side  of  (A20)  and  using  (A22)  for  the  se¬ 
cond  term  on  the  right  hand  side  of  (A20)  we  arrive  at  the  following  expression  for  the  Jacobian: 

u  ,  2a  a  sin^ 0c  - sin^ 0, 

--^  =  /i  +  -tan  0„-cos0^+ - (A23) 


sin  0c 


But  (A23)  can  be  rewritten  as 

h^J  h  2a  -2ACOS0;. +  sin^0c -sin^0.  ..... 

- =  -sec"^  9„  + - ^ ^ ^  .  (A24) 

sin0^  2  "  2h 

However,  the  numerator  of  the  last  term  in  (A24)  vanishes  with  the  result,  after  using  (A17), 
that  the  Jacobian  of  the  transformation  from  ocean  to  sky  is  given  by 


,  _  sin  0^  sec^  9„  sec  O)  sin  0^  sec^  0„ 

J  —  — ■  —  • 

2h  4 


(A25) 
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THE  JACOBIAN  FOR  THE  TRANSFORMATION  BETWEEN  SKY  AND 

SUN 

When  integrating  over  the  solar  disk,  it  is  more  convenient  to  use  a  coordinate  system  re¬ 
ferred  to  the  solar  center  than  a  system  referred  to  the  zenith.  Figure  A-5  shows  the  spherical 
triangle  that  connects  the  two  systems.  From  the  law  of  cosines  and  the  law  of  sines  for  this 
spherical  triangle,  we  have  the  transformation  equations 


cosOs  =  cosSo  cos^  +  sin  sin ^cos;i' 

■  { \  •  e  siny 

sm(0o-<>J  =  sm^— ^  . 

smO, 


(A26) 


The  Jacobian  of  the  transformation  in  (A26)  is  given  by 


J'  = 


dejd^  dejdx 

9<l>sld^  d(f>J9x 


(A27) 


mJ'= 


Calculation  of  the  partial  derivatives  within  the  determinant  gives 

+ cosOq  sin  I  -  sin  cos^cosx  +  sin  6^  sin  sin x\ 

-cos^sin;i<  -sin^cos;!' 

where  m  s  sin205*cos((J)o-  (fe).  When  the  determinant  on  the  right-hand  side  of  (A28)  is  eva¬ 
luated,  we  find  that 

my=(-cos0osin^cos;if  +  sin0ocos^)sin^  . 

From  the  law  of  cosines  for  the  side  0j,  we  have  the  identity 


(A28) 


(A29) 


.  t  COS0, -cos£cos0„ 
sm^cos;i^= - - — ^ — - -  , 


sin^o 

which  can  be  substituted  into  (A29)  to  give 

COS^  -  COS05  cosdf, 


(A30) 


mJ'= 


sin0^ 


sm;i' 


(A31) 


But  the  expression  within  braces  in  (A31)  is  equal  to  sin0j*cos((J)o  -  (fe)  by  the  law  of  co¬ 
sines  for  the  side  |.  Hence,  after  m  is  restored  (A31)  becomes 


sin^  Bs  cos((^o  -  0s sine,  cos(<j>o  -  0jsin;t 


or 


J'- 


smx 

sine^ 


(32) 

(33) 


for  the  Jacobian  of  the  transformation  from  sky  to  sun. 
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UPPER  LIMB 


Z 


Figure  A-5.  Upper  figure:  Coordinates  for  integration  over  the  solar  disk  as  viewed 
looking  up  at  the  sun  from  the  earth.  Lower  figure:  Spherical  triangle  connecting  the 
zenith,  solar  center,  and  a  position  on  the  solar  disk  as  viewed  looking  down  on  the  sun 
toward  the  surface  of  the  earth. 
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THERMAL  EMISSION  FROM  THE  SEA 


In  the  infrared  spectral  region,  solar  effects,  such  as  sun  glint,  become  less  important  and 
thermal  effects,  such  as  black-body  emission,  become  more  important.  The  thermal  emission  of 
the  sea  surface,  which  is  a  major  component  of  the  infrared  radiance  of  the  ocean  horizon,  will 
now  be  derived. 


The  thermal  radiance  emitted  from  a  gray  surface  at  absolute  temperature  T  is 

eN.(T) , 

where  N*{T)  is  Planck  spectral  radiance  of  a  black  body 

..3 


N*(T)=2c^h 


txp{hvl  kT)-\ 


(A34) 

(A35) 


and  E  is  the  emissivity  of  the  object.^  The  energy  arriving  at  sea  water,  an  opaque  object  that 
does  transmit,  can  only  be  reflected  nr  absorbed  according  to  the  law  of  conservation  of 
energy.  From  Kirchoff’s  Law,  the  energy  emitted  is  the  same  as  the  energy  absorbed.  Hence,  for 
sea  water 


e=l-p(ft),v)  .  (A36) 

We  compute  the  thermal  emission  arriving  at  the  receiver  from  the  entire  footprint  shown  in 
figure  9.  The  logical  chain  is  similar  to  that  given  in  section  5  of  the  report;  conservation  of 
radiance  may  not  be  used  directly  and  the  analysis  proceeds  from  the  standpoint  of  power.  The 
thermal  power  emitted  by  a  facet  inside  the  footprint  and  captured  by  the  receiver  is 

dPg-eN*^dB  ,  (A37) 

R 

where  the  quantities  dB,  dS,  and  R  have  previously  been  defined  in  section  5  and  are  also  shown 
in  figure  9.  Integrating  (A37)  over  all  facets  and  using  (7)  and  (A36)  gives 

R 

for  the  total  power  received,  which  can  then  be  divided  by  solid  angle  dS/R^  and  Br,  the  area  of 
the  rough  footprint  projected  toward  the  receiver,  to  give  the  thermal  radiance  emitted  by  the 
sea: 


=N*(Tsea)jl[^-p{(0,v)]qrdS^dSy  .  (A39) 


2  Other  symbols  in  (A3.^)  have  their  usual  meaning. 
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Equation  (A39)  represents  the  thermal  emission  as  an  integral  over  slopes  in  the  ocean.  The 
Jacobian  transforming  from  ocean  to  sky  is  sin  6s»o,  where  a  is  given  by  (24).  Hence, 

nil  In 

=  N*[Tsea )  J  dSs  sin  Os  J  -  fy{(0,  v)]q,  a  (A40) 

0  0 

is  an  alternate  representation  for  the  sea  emission  in  which  integration  occurs  over  the  sky. 
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Cox,  C.,  and  W.  Munk.  1955.  “Some  Problems  in  Optical  Octanography,"  Journal  of  Marine 
Research,  vol.  14,  p.  74. 


A-13 


TR  1660 


APPENDIX  B 

FORTRAN  SOURCE  CODE  LISTINGS 


B-l 


The  following  FORTRAN  programs  are  based  on  the  equations 
presented  in  this  report.  The  equation  being  computed  is  noted  in  columns  73 
through  80  of  each  listing.  Each  program  assembles  without  errors  or 
warnings  on  a  Lahey  F77L/EM-32  compiler  version  5.01  and  runs  on  a 
personal  computer  provided  with  a  Lahey/Phar  Lap  386|DOS-Extender. 

These  programs  accept  input  from  the  keyboard  and  produce  output 
either  to  the  screen  or  to  the  hard  drive.  In  the  latter  case  output  consists  of  a 
file  in  ASCn  format  which  can  be  read  by  commercial  software  plotting 
packages.  An  example  of  file  output  follows  those  listings  which  created 
them.  Many  of  the  figures  in  this  report  were  generated  with  these  programs 
and  GRAFTOOL™. 

For  a  copy  of  this  source  code  on  3  1/2  inch  disk  please  write  to  the 
author. 
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SOURCE  CODE 
TABLE  OF  CONTENTS 


Function 

Rho 

Sky 

Sun 

Tilt 


Ellipse 

Plotglnt 

Plotp 

Plotq 

Plotrcvr 

Plotsun 

Plotwind 


Testb 

Testp 

Testq 

Testsea 

TestslQ^ 

Testsun 

Testilt 


Functions  and  Subroutines 


Cox-Munk  functions. 

Reflectivity  of  sea  water. 

Reflected  sky  radiance  with  sea  emission. 

Glint  radiance  ratio. 

Facet  properties. 

Plot  Drivers 

ASCn  file  for  edge  of  tolerance  ellipse. 

ASCn  file  for  glint  versus  pixel  position. 

ASCn  file  for  p  versus  slope. 

ASCn  file  for  q  versus  slope. 

ASCn  file  for  glint  versus  receiver  zenith  angle. 
ASCn  file  for  glint  versus  solar  zenith  angle. 
ASCn  file  for  glint  versus  wind  speed. 


Checks  function  BoverA  firom  keyboard. 
Checks  function  p  fi-cnn  keyboard. 

Checks  function  q  fi'om  keyboard. 
CoX'Munk  keyboard  and  input  file  driver. 
Checks  "Sky"  fi’om  keyboard. 

Checks  "Sun"  fi-om  keyboard. 

Checks  "Tilt"  fi’om  keyboard. 


oooo  oooooo  o  o  n  oo  o  o  oo  o'o  no 


C:\COX\FUNCTION.FOR  8/1/94 


C 

*****************************  FUNCTION. FOR  ***************************** 

*  * 

*  Various  useful  Cox-Munk  functions.  * 

*  * 

*  Latest  revision;  August  1,  1994  * 

*  * 
************************************************************************ 

C 

FUNCTION  Su(W) 

gives  the  rms  upwind  slope  deviation: 

Su  =  SQRT(3.16E-3*W)  (2) 

END 

FUNCTION  SC(W) 

gives  the  rms  crosswind  slope  deviation: 

Sc  =  SQRT(3.E~3  +  1.92E-3*W)  (2) 

END 

FUNCTION  p(Sx,  Sy,  W) 

is  the  Cox-Munk  wave  slope  probability  distribution. 

COMMON/Constants/Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

pO  =  l./(2.*Pi*Su(W)*Sc(W)) 

Arg  =  ( (Sx/Su(W) )**2  +  (Sy/Sc(W) )**2)/2. 

IF  ((Arg)  .GE.  LOG(pO/Delta) )  THEN 
p  =  Zero 

ELSE 

p  =  pO*EXP(-Arg)  (2) 

END  IF 
END 

FUNCTION  BoverA(Tx,  Px,  W) 

is  area  B  of  the  ergodic  cap  projected  toward  u  =  (Tx,Px) 
divided  by  the  area  A  of  the  cap  projected  vertically. 

With  N  =  2,  M  =  7:  Accuracy  is  <  1%  for  W  below  40  m/s  and 
>  1%  for  W  above  40  m/s  (due  to  failure  to  abandon  Cox-Munk 
approximation  for  the  integral  above  40  m/s) . 

Times  for  a  486  personal  computer  with  a  50  MHz  clock: 

Cox-Munk  region,  60  microseconds;  integral  region,  5  milli¬ 
seconds  . 

COMMON/Constants/Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

N  =2 
M  =7 
C 

Sav  =  (Su(W)  +  Sc(W))/2. 
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Sntax  =  N*2.303*Sav 

ds  =  Smax/M 

C 

A  =  SIN(TX) *C0S(PX)  (A2) 

B  =  SIN(TX) *SIN(Px) 

C  =  COS(TX) 

C 

IF  (Tx  .LE.  D2R*58.)  THEN 


BoverA  = 

C 

(A9) 

ELSE  IF  ( (Tx 
BoverA  = 

.LE.  D2R*64.) 

C 

.AND. 

(W  .LE. 

26.)) 

THEN 

(A9) 

ELSE  IF  ((Tx 
BoverA  = 

.LE.  D2R*68.) 
C 

.AND. 

(W  .LE. 

10.)) 

THEN 

(A9) 

ELSE  IF  ((Tx 
BoverA  = 

.LE.  D2R*72.) 
C 

.AND. 

(W  .LE. 

11.)) 

THEN 

(A9) 

ELSE  IF  ((Tx 
BoverA  = 

.LE.  D2R*76.) 

C 

.AND. 

(W  .LE. 

6.)) 

THEN 

(A9) 

ELSE  IF  ( (Tx 
BoverA  = 

.LE.  D2R*80.) 
C 

.AND. 

(W  .LE. 

2.)) 

THEN 

(A9) 

ELSE 

SUM  =  0. 

DO  Sx  =  -Smax,  Smax,  dS 

Symax  =  SQRT(ABS(Smax**2  -  Sx**2)) 

DO  Sy  =  -Symax,  Symax,  dS 

Ratio  =  -  A*Sx  -  B*Sy  +  C  (A8) 

IF  (Ratio  .GT.  0.)  SUM  =  Ratio*p(Sx,  Sy,  W)  +  SUM  (A7) 
END  DO 
END  DO 

BoverA  =  SUM*dS**2 

END  IF 
C 

END 

C 

FUNCTION  qo(Tx,  Px,  Sx,  Sy,  W) 

C  is  the  interaction  probability  density  without  "BoverA”. 

C 

A  =  SIN(Tx)*COS(Px)  (A2) 

B  =  SIN(TX)*SIN(PX) 

C  =  COS(TX) 

C 

Ratio  =  -  A*Sx  -  B*Sy  +  C  (A8) 

IF  (Ratio  .GT.  0.)  THEN 

qo  =  Ratio*p(Sx,  Sy,  W)  (13) 

ELSE 

qo  =  0.  (13) 

END  IF 
C 

END 

C 

FtJNCTION  q(Tx,  Px,  Sx,  Sy,  W) 

C  is  the  interaction  probability  density. 

C 

q  =  qo(Tx,  Px,  Sx,  Sy,  W) /BoverA(Tx,  Px,  W)  (13) 

C 

END 

C 

FUNCTION  Ns(a,b,c,T) 


Page  2 


oooo  o  o  o  ooo 


c 


c 

c 

c 
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is  an  analytic  form  for  sky  radiance  (W  m-2  sr-1  (cm-l)-i) 
as  a  function  of  zenith  angle  T  (radians) . 

REAL* 4  Ns 

Ns  =  a  +  b*exp(T/c) 

END 

FUNCTION  BB(V,  T) 

is  the  spectral  radiance  of  a  black  body  in  units  of 
W  m-2  sr-1  (cm-l)-l.  V  is  the  wave  number  in  cm-1  and 
T  is  the  temperature  in  Kelvin. 

Cl  =  11910.62 
C2  =  14387.86 

W  =  1.E4/V 
A  =  Cl/W**3 
X  =  C^/(W*T) 

BB  =  A/(Exp(X)  -  1.) 

RETURN 

END 
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FUNCTION  Rho( Omega,  V) 

C 

COMMON  /Sealndex/  AlphaOl (100) ,  Alpha02(20), 

+  BetaOl  (100),  Beta02  (20) 

C 

ititititidtit'kifk'k’kifk'kit'k'k’k'k'kie'k'k’k'kifkifk  RHO •  FOR 

it  it 

*  Calculates  reflectivity  of  sea  water  between  52.63  cm-1  and  * 

*  25,000  cm-1  using  equations  (74)  and  (78)  of  Stratton,  "Electro-  * 

*  magnetic  Theory",  1941,  page  505,  ff.  The  sea  water  is  assumed  to  * 

*  be  a  conducting  medium;  both  real  and  imaginary  parts  of  the  * 

*  index  of  water  are  used.  The  notation  of  Stratton  is  adhered  to  * 


*  as  closely  as  possible.  * 

*  * 

*  Omega  is  the  angle  of  incidence  in  radians;  V  is  the  wave-  * 

*  number  in  cm-1;  Rho  is  the  reflectivity.  * 

*  * 

*  Last  revision:  March  10,  1994  * 

*  * 


c 

C  The  four-point  interpolation  functions  are: 

WM(X)  =  (X  -  l.)*(X  -  2.)*X/6. 

W0(X)  =  (X  -  1,)*(X  -  2.)*(X  +  l.)/2. 

W1(X)  =  (X  +  l.)*(X  -  2.)*X/2. 

W2(X)  =  (X  +  l.)*(X  -  l.)*X/6. 

C 

IF  (V  .EQ.  0.)  THEN 

C  set  reflectivity  to  100%  and  return: 

Rho  =  1. 

RETURN 
END  IF 
C 

W  =  1.E4/V 
C 

IF  (W  .LT.  0.399999)  THEN 
C  print  error  messate  and  return: 

Rho  =  0. 

WRITE  (6,  1000)  V 
RETURN 
END  IF 
C 

IF  (0.4  .LE.  W  .AND.  W  .LE.  19.8)  THEN 
C  use  Lagrange  4  point  interpolation  on  Block  01  data  which 

C  are  at  0.2  urn  spacing  between  0.2  and  20  urn: 

I  =  INT(W/0.2) 

Fr  =  MOD(W,  0.2)/0.2 

Alphal  =  W2 (Fr) * AlphaOl (I  +  2)  -  Wl(Fr) *Alpha01 (I  +  1) 

+  +  W0(Fr)*Alpha01(I)  -  WM(Fr) *Alpha01 (I  -  1) 

Betal  =  W2(Fr)*Beta01  (I  +  2)  -  Wl (Fr) *Beta01  (I  +  1) 

+  +  W0(Fr)*Beta01  (I)  -  WM(Fr) *Beta01  (I  -  1) 

END  IF 
C 

IF  (19.8  .LT.  W  .AND.  W  .LT.  190.)  THEN 
C  use  Lagrange  4  point  interpolation  on  Block  02  data  which 

C  are  at  10  urn  spacing  between  10  and  200  urn: 

I  =  INT(W/10.) 
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C 

C 


C 


Fr  =  MOD(W,  10.)/10. 

Alphal  =  W2(Fr)*Alpha02(I  +  2) 

+  +  W0(Fr)*Alpha02(I) 

Betal  =  W2(Fr)*Beta02  (I  +  2) 

+  +  W0(Fr)*Beta02  (I) 

END  IF 


Wl(Fr)*Alpha02 (I  +  1) 
WM(Fr)*Alpha02(I  -  1) 
Wl(Fr)*Beta02  (I  +  1) 
WM(Fr)*Beta02  (I  -  1) 


IF  (190.  .LE.  W)  THEN 
print  error  message  and  return: 
RHO  =  0. 

WRITE  (6,  1000)  V 
RETURN 
END  IF 


G  *  Alphal**2  -  Betal**2  -  SIN(Omega) **2 
H  *  4*{Alphal**2) *(Betal**2) 

P  =  SQRT(0.5*(-G  +  SQRT(H  +  G**2))) 

Q  =  SQRT(0.5*(+G  +  SQRT(H  +  G**2))) 

Stratton,  Equation  (74),  p.  505: 

C  =  (Q  -  COS (Omega) ) **2  +  P**2 
D  =  (Q  +  COS (Omega) ) **2  +  P**2 
Rp  =  C/D 

Stratton,  Equation  (77),  p.  506: 

E  =  ( (Alphal<r*2  -  Betal**2)  *COS  (Omega)  -  Q)  **2 

F  =  ( (Alphal**2  -  Betal**2) *COS (Omega)  +  Q) **2 

T  =  (2*Alphal*Betal*COS(Omega)  -  P)**2 

U  =  (2*Alphal*Betal*COS (Omega)  +  P)**2 

Rs  =  (E  +  T)/(F  +  U) 

Rho  =  (Rp  +  Rs)/2. 

RETURN 


1000  FORMAT  ('  *****  WARNING  -  INPUT  FREQUENCY  =  ',  1PG12.5,  'CM-1', 

+  /,  '  OUTSIDE  VALID  RANGE  OF  52.63  TO  25,000  CM-1  ******',/) 


END 


BLOCK  DATA  SEA 
C 

COMMON  /Sealndex/  AlphaOl (100) ,  Alpha02(20), 

+  BetaOl  (100) ,  Beta02  (20) 

C 

*  * 

*  These  data  for  the  optical  index  of  water  have  been  taken  from  * 

*  the  literature.  From  0.2  to  2  microns  (blocks  01  up  to  second  * 

*  entry  of  row  B)  and  from  10  to  200  microns  (blocks  02)  the  data  * 

*  are  from  G.  M.  Hale  and  M.  R.  Querry,  "Optical  Constants  of  Water* 

*  in  the  200-nm  to  200-um  Wavelength  Region,"  Appl.  Opt.  3,  555  * 

*  (1973) .  These  data  are  for  pure  water.  * 
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*  * 

*  From  2.2  to  20  microns  (blocks  01  from  the  third  entry  of  row  B  * 

*  to  the  end)  the  data  are  from  M.  R.  Querry,  W.  E.  Holland,  R.  C.  * 

*  Waring,  L.  M.  Earls,  and  M.  D.  Querry,  "Relative  Reflectance  and  * 

*  Complex  Refractive  Index  in  the  Infrared  for  Saline  Environmental* 

*  Waters,"  J,  Geophys.  Res.  82,  1425  (1977),  Table  3,  Pacific  * 

*  Ocean  colvunns.  These  data  are  for  salt  water.  * 

*  * 
************************************************************************ 
C 

.C  Real  part  of  the  index  of  sea  water  from  0.2  to  20  microns  in 
C  steps  of  0.2  microns: 

DATA  AlphaOl  / 


A 

1.396, 

1.339, 

1.332, 

1.329, 

1.327, 

1.324, 

1.321, 

1.317 

B 

1.312, 

1.306, 

1.303, 

1.287, 

1.251, 

1.151, 

1.384, 

1.479 

C 

1.421, 

1.388, 

1.368, 

1.355, 

1.347, 

1.339, 

1.335, 

1.335 

D 

1.332, 

1.324, 

1.312, 

1.296, 

1.268, 

1.271, 

1.371, 

1.353 

E 

1.340, 

1.330, 

1.324, 

1.319, 

1.314, 

1.307, 

1.302, 

1.297 

F 

1.291, 

1.286, 

1.279, 

1.272, 

1.268, 

1.264, 

1.258, 

1.249 

G 

1.240, 

1.229, 

1.218, 

1.204, 

1.190, 

1.177, 

1.165, 

1.151 

H 

1.140, 

1.132, 

1.124, 

1.119, 

1.121, 

1.126, 

1.134, 

1.142 

I 

1.152, 

1.164, 

1.177, 

1.189, 

1.201, 

1.212, 

1.224, 

1.234 

J 

1.242, 

1.253, 

1.261, 

1.273, 

1.284, 

1.296, 

1.309, 

1.320 

K 

1.331, 

1.339, 

1.349, 

1.358, 

1.366, 

1.379, 

1.390, 

1.399 

L 

1.408, 

1.417, 

1.426, 

1.435, 

1.443, 

1.450, 

1.458, 

1.464 

M 

1.470, 

1.474, 

1.477, 

1.480 

/ 

C 

C  Real  part  of  the  index  of  sea  water  from  10  to  200  microns  in 
C  steps  of  10  microns: 

DATA  Alpha02  / 

A  1.218,  1.480,  1.551,  1.519,  1.587,  1.703,  1.821,  1.886, 

B  1.924,  1.957,  1.966,  2.004,  2.036,  2.056,  2.069,  2.081, 

C  2.094,  2.107,  2.119,  2.130  / 

C 

C  Imaginary  part  of  the  index  of  sea  water  from  0.2  to  20  microns  in 
C  steps  of  0.2  microns: 


DATA 

BetaOl 

/ 

A 

0.000, 

0.000, 

0.000, 

0.000, 

0.000, 

0.000, 

0.000, 

0.000, 

B 

0.000, 

0.001, 

0.000, 

0.001, 

0.003, 

0.114, 

0.263, 

0.085, 

C 

0.018, 

0.005, 

0.003, 

0.004, 

0.007, 

0.011, 

0.016, 

0.016, 

D 

0.013, 

0.011, 

0.010, 

0.013, 

0.032, 

0.108, 

0.087, 

0.044, 

E 

0.035, 

0.033, 

0.032, 

0.031, 

0.031, 

0.032, 

0.032, 

0.033, 

F 

0.034, 

0.036, 

0.038, 

0.041, 

0.044, 

0.046, 

0.046, 

0.047, 

G 

0.048, 

0.050, 

0.054, 

0.060, 

0.068, 

0.079, 

0.091, 

0.107, 

H 

0.125, 

0.145, 

0.166, 

0.191, 

0.216, 

0.239, 

0.260, 

0.279, 

I 

0.297, 

0.313, 

0.326, 

0.338, 

0.347, 

0.357, 

0.363, 

0.371, 

J 

0.377, 

0.385, 

0.393, 

0.401, 

0.407, 

0.413, 

0.417, 

0.418, 

K 

0.420, 

0.422, 

0.424, 

0.427, 

0.430, 

0.432, 

0.432, 

0.432, 

L 

0.431, 

0.430, 

0.429, 

0.427, 

0.425, 

0.423, 

0.420, 

0.416, 

M 

0.414, 

0.410, 

0.406, 

0.393 

/ 

C 

C  Imaginary  part  of  the  index  of  sea  water  from  10  to  200  microns  in 
C  steps  of  10  microns: 

DATA  Beta02  / 

A  0.051,  0.393,  0.328,  0.385,  0.514,  0.587,  0.576,  0.547, 

B  0.536,  0.532,  0.531,  0.526,  0.514,  0.500,  0.495,  0.496, 

C  0.497,  0.499,  0.501,  0.504  / 
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C 


END 
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SUBROUTINE  Sky ( V , W , Nsky , Nsea ) 

C 

COMMON/ IN/  Ts,  Ps,  Tr,  Pr 
COMMON/OUT/  Sx,  Sy,  Omega,  Tn,  Pn 
COMMON/ Constants/  Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

COMMON/Sea/  a,  b,  c,  Tsea 
C 


*  * 

*  Returns  reflected  sky  radiance  Nsky  by  integrating  the  source  * 

*  radiance  Ns  over  the  sky  dome.  Returns  emitted  sea  radiance  * 

*  Nsea  by  integrating  the  sea  emission  BB  over  the  sky  dome.  * 

*  * 

*  V  is  the  wavenumber  at  which  reflection  occurs  and  W  is  the  * 

*  wind  speed  at  the  time  of  reflection  and  emission.  * 

*  * 

*  USES  Function,  Tilt,  Rho  * 

*  * 

*  Latest  revision;  August  1,  1994  * 

*  * 


************************************************************************ 

c 

REAL*4  Ns,  Nsky,  Nsea 
C 

N  =  18 

C  for  5  degree  sky  segments  and  1%  precision. 

C 

SumSky  *  0. 

SumSea  =  0. 

dTs  =  (Pi/2.) /FLOAT (N) 

DO  I  =  1,  N 

Ts  =  (I  -  0.5)*dTs 

K  =  INT(2.*Pi*Sin(Ts)/dTs  +  0.4999) 

dPs  =  2*Pi/FLOAT(K) 

SumPsSky  =  0. 

SumPsSea  =  0. 

DO  J  =  1,  K 

Ps  =  -  Pi  +  (J  -  0.5)*dPs 
CALL  Tilt 
SumPsSky  = 

+  dPs*Rho ( Omega, V) /COS(Tn) **4*p(Sx,Sy,W) 

SumPsSea  = 

+  dPs*(l  -  Rho(Omega,V))/COS(Tn)**4*p(Sx,Sy,W) 

END  DO 

SumSky  =  Sin(Ts) *Ns(a,b,c,Ts)*SumPsSky  +  SumSky 
SumSea  =  Sin(Ts) *BB(V,  Tsea) *SumPsSea  +  SumSea 
END  DO 

Nsky  =  SumSky*dTs/4./BoverA(Tr,Pr,W)  (30) 

Nsea  =  SumSea*dTs/4 . /BoverA(Tr,Pr,W)  (A40) 

C 

RETURN 

C 

END 


+  SumPsSky  (30) 

+  SumPsSea  (A40) 

(30) 

(A40) 
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SUBROUTINE  Sun(V,  W,  NoverNo) 

C 

COMMON/ IN/  Ts,  Ps,  Tr,  Pr 
COMMON/OUT/  Sx,  Sy,  Omega,  Tn,  Pn 
COMMON/ Constants/  Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

COMMON/ Sea/  a,  b,  c,  Tsea 
C 

*********************************  Sun • FOR  ****************************** 
*  * 

*  Exact  calculation  of  sun  glint  ratio  N/No.  Uses  rectangular  * 

*  coordinates  (x,  y)  over  disk  rather  than  radial  coordinates  * 

*  (Xi,  Chi)  given  in  report.  * 

*  * 

*  V  is  the  wavenumber  in  cm-1,  W  is  the  wind  speed  in  m  s-1,  and  * 

*  NoverNo  is  the  glint  ratio,  an  effective  reflectivity.  This  * 

*  subroutine  should  be  called  with  the  common  block  "IN”  loaded  * 

*  with  Ts,  Ps  the  location  of  the  solar  center  and  Tr,  Pr  the  * 

*  location  of  the  receiver.  * 

*  * 

*  The  larger  the  value  of  M,  the  y  coordinate  step  size,  the  more  * 

*  precise  and  slower  the  sum.  For  fixed  M  precision  improves  * 

*  with  wind  speed.  For  W  =  l  m  s-1  and  M  =  5,  the  precision  is  * 

*  better  than  1.5  %  around  the  center  of  the  glint  pattern  until  * 

*  the  receiver  zenith  angle  exceeds  89.5  degrees.  * 

it  it 

*  USES  Tilt,  Function,  Rho  * 

*  Latest  revision:  August  12,  1994  * 

*  Bug:  Divides  by  zero  when  To  =  0  (sun  is  on  the  zenith) .  * 

*  * 
************************************************************************ 
C 

To  =  Ts 
Po  =  Ps 
C 

REAL* 4  NoverNo 
M  =  5 

N  =  2*M  +  1 


SUM  =  0. 

dY  =  2.*Epsilon/N 
DO  I  =  1,  N 

Y  =  Epsilon  -  (I  -  0.5) *dY 
Xmax  =  SQRT(Epsilon**2  -  Y**2) 

K  =  INT(2.*Xmax/dY  +  0nem*0.5) 

dX  =  2.*Xmax/(K) 

DO  J  =  1,  K 

X  =  Xmax  -  (J  -  0.5) *dX 

Ts  =  To  -  Y 

Ps  =  Po  -  X/SIN(To) 

CALL  Tilt 

SUM  =  Rho(Omega,  V) / (COS(Tn) **4) *p(Sx,  Sy,  W) *dX  +  SUM 
END  DO 
END  DO 

NoverNo  =  SUM*dY/BoverA(Tr,  Pr,  W)/4. 

Ts  =  To 
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C;\COX\SUN.FOR  8/12/94 


C 

C 


Ps  »  Po 

RETURN 

END 


Page  2 


\COX\TILT.FOR  8/1/94 


SUBROUTINE  TILT 


COMMON/ IN/Ts , Ps , Tr , Pr 

COMMON/OUT/ Sx , Sy , Onega , Tn , Pn 

COMMON/ Constants/ Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

C 

*********************************  TILT. FOR  ***************************** 
*  * 

*  Computes  various  features  of  facet  geometry  given  the  * 

•  *  source  and  receiver  coordinates.  * 

*  * 

*  Angles  passed  in  radians.  Azimuths  +  CCW  looking  along  nadir.  * 

It  it 

*  Latest  revision:  August  1,  1994  * 

*  * 

************************************************************************ 

c 

IF  ((Ts  .GT.  Pi/2.*Onep)  .OR.  (Tn  .GT.  Pi/2.*Onep))  THEN 

PRINT  *,  "Error  from  'TILT’:  Zenith  angle  of  source  or" 
PRINT  *,  "  receiver  exceeds  90  degrees." 

RETURN 

END  IF 


As  =  SIN(Ts)*COS(Ps) 

Ar  =  SIN(Tr)*COS(Pr) 

Bs  =  SIN(Ts)*SIN(Ps) 

Br  =  SIN(Tr)*SIN(Pr) 

Cs  *  COS(Ts) 

Cr  =  COS(Tr) 

IF  (ABS(As  +  Ar)  .LE.  Delta)  THEN 
Sx  =  0. 

ELSE  IF  ((Cs  +  Cr)  ,LE.  Delta)  THEN 
Sx  =  SIGN (Infinity,  -(As  +  Ar) ) 

ELSE 

Sx  =  -  (As  +  Ar)/(Cs  +  Cr) 

END  IF 


(A12) 


IF  (ABS(Bs  +  Br)  .LE.  Delta)  THEN 
Sy  =  0. 

ELSE  IF  ((Cs  +  Cr)  .LE.  Delta)  THEN 
Sy  =  SIGN (Infinity,  -(Bs  +  Br) ) 

ELSE 

Sy  =  -  (Bs  +  Br)/(Cs  +  Cr) 

END  IF 


(A13) 


D  =  (1.  +  As*Ar  +  Bs*Br  +  Cs*Cr)/2. 

IF  (D  .LE.  Delta)  D  =  0. 

S  =  SQRT  (D) 

If  ((Onem  .LE.  S)  .AND.  (S  .LE.  Onep))  THEN 
Omega  =  0. 

ELSE  IF  ((-Onep  .LE.  S)  .AND.  (S  .LE.  -Onem))  THEN 
Omega  =  Pi 

ELSE 

Omega  =  ACOS(S) 

END  IF 


(A14) 

(A14) 


(A14) 
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C; \COX\TILT.FOR  8/1/94 


C 

Tn  »  ATAN(SQRT(SX**2  +  Sy**2)) 

C 

IF  ((ABS(Sx)  .LE.  Delta)  .AND.  (ABS(Sy)  .LE.  Delta))  THEN 
Pn  ®  0. 

ELSE  IF  ((Cs  +  Cr)  .LE,  Delta)  THEN 
Pn  «  (Ps  +  Pr)/2. 

ELSE 

Pn  *=  ATAN2(-Sy,  -Sx) 

END  IF 
C 

RETURN 

C 

END 


(A6) 


!  (A2)  +  (A5) 
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C:\COX\ELLIPSE.FOR  8/1/94 


COMMON/ IN/Ts , Ps , Tr , Pr 
COMMON/OUT/Sx , Sy , Omega , Tn, Pn 
COMMON/Constants/Pl,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

C 

****************************  ELLIPSE. FOR  ******************************* 


*  Calculates  the  Cox-Munk  ellipse  in  slope  space  such  that  the  * 

*  slopes  within  the  ellipse  will  reflect  the  sun's  disk  into  the  * 

*  receiver.  Uses  spherical  triangle  defined  by  zenith,  solar  * 

■*  center,  and  edge  of  solar  disk  to  compute  slopes  as  the  para-  * 

*  meter  "Xi"  is  stepped  around  the  edge  of  the  solar  disk.  The  * 

*  spherical  triangle  is  shown  in  figure  A-5.  * 

*  * 

*  Computes  value  of  q,  the  interaction  prob.  density,  for  edge  * 

*  of  the  ellipse  and  prints  it  to  a  file  in  Sx,  Sy,  q  form  for  * 

*  plotting  in  three  dimensions.  These  data  define  the  edge  of  * 

*  the  glint  cap  which  is  the  top  of  a  glint  column.  * 

*  * 

*  USES  TILT,  FUNCTION  * 

*  Last  revised:  August  1,  1994  * 

*  * 
************************************************************************ 


C 

CHARACTER*80  H 
H(l:80)  =  "  " 

H(9:ll)  =  "\SX" 

H(23:24)  =  "Sy" 

H(36:36)  =  "q" 

H(44:51)  =  "Sx  -  Sxo" 

H(58;65)  =  "Sy  -  Syo" 

H(75:76)  =  "Xi" 

C 

Pi  =  4.*ATAN(1.) 

R2D  =  180. /Pi 

D2R  =  Pi/180. 

Epsilon  =  D2R*0.2659 
Delta  =  1.4E-6 
Onem  =  1.  -  Delta 
Onep  =  1.  +  Delta 
Infinity  =  999999 
Zero  =  0. 

TO  =  273.15 
C 

IPR  =  10 

OPEN  (IPR,  CARRIAGE  CONTROL  =  'FORTRAN',  FILE  =  'E.DAT') 
C 

PRINT  *,  "CURRENT  WIND  SPEED  (m  s-1)  =  " 

READ  (5,  *)  W 
IF  (W  .LT.  0.01)  THEN 
W  =  0.01 
END  IF 

PRINT  *,  "CAMERA  ZENITH  ANGLE,  AZIMUTH  (deg)  =  " 

READ  (5,  *)  Tr,  Pr 

PRINT  *,  "SOLAR  ZENITH  ANGLE,  AZIMUTH  (deg)  =  " 

READ  (5,  *)  Tso,  Pso 
C 
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C:\COX\ELLIPSE.FOR  8/1/94 


Tr  »  D2R*Tr 
Pr  »  D2R*Pr 
TSO  »  D2R*Tso 
Pso  »  D2R*Pso 
C 

Ts  =  Tso 
Ps  =  Pso 
CALL  TILT 
Sxo  -  Sx 
Syo  =  Sy 

Area=  Pi* (Epsilon) **2/ (4 . *COS (Omega) *COS(Tn) **3)  !(24)+(25) 

C 

WRITE  (IPR,  90)  "\TOLERANCE  CAP:  WIND  SPEED  (m/s)  =  •' ,  W 
WRITE  (IPR,  80)  "XCAMERA  ZEN,  AZM  (deg)  =  ",  R2D*Tr,  R2D*Pr 
WRITE  (IPR,  80)  "\SOLAR  ZEN,  AZM  (deg)  =  ",  R2D*Tso,  R2D*Pso 
WRITE  (IPR,  90)  "\CENTRAL  INCIDENCE  ANGLE  (deg)  =  ",  R2D*Omega 
WRITE  (IPR,  90)  "\CENTRAL  TILT  (deg)  =  ",  R2D*Tn 
WRITE  (IPR,  70)  "\ APPROXIMATE  AREA  =  ",  Area 
WRITE  (IPR,  60)  "\CENTRAL  SLOPES  Sxo,  Syo  =  ",  Sxo,  Syo 
WRITE  (IPR,  20)  H 
C 

IF  ((Tso  .GE.  Onem*Epsilon)  .AND.  (Tso  .LE.  Onep*Epsilon) ) 

C  reset  zenith  angle  to  value  in  stable  range: 

+  Tso  =  Tso*0nep 
A  =  COS (Tso) *COS (Epsilon) 

B  «=  SIN  (Tso)  *SIN  (Epsilon) 

OldC  =  -  Delta 

DO  Xi  =  D2R*0.,  D2R*360.,  D2R*10. 

Ts  =  ACOS(A  +  B*C0S(Xi))  (A26) 

C  =  SIN (Epsilon) /SIN (Ts)*SIN(Xi)  (A26) 

IF  (C  .GE.  +  Onem)  C  =  +  Onem 
IF  (C  .LE.  -  Onem)  c  =  -  Onem 
Ps  =  Pso  +  ASIN(C) 

IF  (Tso  .LT.  Onem* Epsilon)  THEN 

C  zenith  pierces  solar  disk  and  azimuth  cycles  +  pi  to  -  Pi: 

IF  ((C  .GE.  0.)  .AND.  (C  .GE.  OldC))  THEN 
Ps  =  -  Ps  +  Pi 

ELSE  IF  ((C  .LT.  0.)  .AND.  (C  .GE.  OldC))  THEN 
Ps  =  -  Ps  -  Pi 

END  IF 
END  IF 
OldC  =  C 
CALL  TILT 

WRITE  (IPR,  50)  Sx,  Sy,  q(Ts,  Ps,  Sx,  Sy,  W) , 

+  Sx  -  Sxo,  Sy  -  Syo,  R2D*Xi 

END  DO 
C 

PRINT  *,  "DONE.  DATA  IN  'E.DAT'" 

C 

20  FORMAT  (A80,  /) 

50  FORMAT  (2(1PE14.3),  0PF10.3,  2(1PE14.3),  OPFIO.O) 

60  FORMAT  (A40,  2(1PE15.3),  /) 

70  FORMAT  (A40,  1PE15.3,  /) 

80  FORMAT  (A40,  2F9.2,  /) 

90  FORMAT  (A40,  F9 . 2 ,  /) 

C 

END 
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C:\COX\E.DAT  4/28/94 


\TOLERANCE  CAP;  WIND  SPEED  (m/s)  = 
\CAMERA  ZEN,  AZM  (deg)  = 
\ SOLAR  ZEN,  AZM  (deg)  = 
\ CENTRAL  INCIDENCE  ANGLE  (deg)  = 
\ CENTRAL  TILT  (deg)  = 
\APPROXIMATE  AREA  = 


\ CENTRAL  SLOPES  Sxo,  Syo  = 


\SX 

Sy 

q 

O.OOOE+00 

2.320E-03 

5.666 

2.291E-03 

2.286E-03 

5.667 

4.515E-03 

2.184E-03 

5.668 

6.606E-03 

2.017E-03 

5.669 

8.504E-03 

1.790E-03 

5.672 

1.015E-02 

1.509E-03 

5.675 

1.150E-02 

1.183E-03 

5.680 

1.250E-02 

8.205E-04 

5.686 

1.313E-02 

4.325E-04 

5.692 

1.336E~02 

3.055E-05 

5.700 

1.319E-02 

-3.733E-04 

5.709 

1.261E-02 

-7.667E-04 

5.719 

1.165E-02 

-1.137E-03 

5.729 

1.032E-02 

-1.474E-03 

5.738 

8.677E-03 

“1.765E-03 

5.747 

6.758E-03 

-2.002E-03 

5.755 

4.628E-03 

-2.177E-03 

5.760 

2.351E-03 

-2.284E-03 

5.764 

O.OOOE  00 

-2.321E-03 

5.765 

-2.351E-03 

-2.284E-03 

5.764 

-4.627E-03 

-2.177E-03 

5.760 

-6.758E-03 

-2.002E-03 

5.755 

“8.677E-03 

-1.765E-03 

5.747 

-1.032E-02 

-1.474E-03 

5.738 

-1.165E-02 

-1.137E-03 

5.729 

-1.261E-02 

-7.667E-04 

5.719 

-1.319E-02 

-3.733E-04 

5.709 

-1.336E-02 

3.055E-05 

5.700 

-1.313E-02 

4.325E-04 

5.692 

-1.250E-02 

8.205E-04 

5.686 

-1.150E-02 

1.183E-03 

5.680 

-1.015E-02 

1.509E-03 

5.675 

-8.504E-03 

1.790E-03 

5.672 

-6.606E-03 

2.017E-03 

5.669 

-4.514E-03 

2.184E-03 

5.668 

-2.291E-03 

2.286E-03 

5.667 

O.OOOE+00 

2.320E-03 

5.666 

1C  .  00 

80.00  270.00 

80.00  90.00 

80.00 
0.00 

9.741E-05 


C-uOOE+00  O.OOOE+00 


Sx  -  Sxo 

Sy  -  Syo 

Xi 

O.OOOE+00 

2.320E-03 

0. 

2.291E-03 

2.286E-03 

10. 

4.515E-03 

2. 184E-03 

20. 

6.606E-03 

2.017E-03 

30. 

8.504E-03 

1.790E-03 

40. 

1.015E-02 

1.509E-03 

50. 

1.150E-02 

1.183E-03 

60. 

1.250E-02 

8.205E-04 

70. 

1.313E-02 

4.325E-04 

80. 

1.336E-02 

3.055E-05 

90. 

1.319E-02 

-3.733E-04 

100. 

1.261E-02 

-7.667E-04 

110. 

1.165E-02 

-1.137E-03 

120. 

1.032E-02 

“1.474E-03 

130. 

8.677E-03 

-1.765E-03 

140. 

6.758E-03 

-2.002E-03 

150. 

4.628E-03 

-2.177E-03 

160. 

2.351E-03 

-2.284E-03 

170. 

O.OOOE+00 

-2.321E-03 

180. 

-2.351E-03 

-2.284E-03 

190. 

-4.627E-03 

-2.177E-03 

200. 

-6.758E-03 

-2.002E-03 

210. 

-8.677E-03 

-1.765E-03 

220. 

-1.032E-02 

-1.474E-03 

230. 

-1.165E-02 

-1.137E-03 

240. 

-1.261E-02 

“7.667E-04 

250. 

-1.319E-02 

-3.733E-04 

260. 

-1.336E-02 

3.055E-05 

270. 

-1.313E-02 

4.325E-04 

280. 

-1.250E-02 

8.205E-04 

290. 

-1.150E-02 

1.183E-03 

300. 

-1.015E-02 

1.509E-03 

310. 

-8.504E-03 

1.790E-03 

320. 

-6.606E-03 

2.017E-03 

330. 

-4.514E-03 

2.184E-03 

340. 

-2.291E-03 

2.286E-03 

350. 

O.OOOE+00 

2.320E-03 

360. 
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C:\COX\PLOTGLNT.FOR  4/28/94 


C 

COMMON/ IN/To , Po , Tr , Pr 
COMMON/OUT/ Sx , Sy , Omega , Tn , Pn 
COMMON/Constants/Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

C 

******************************  PlotGlnt.FOR  **************************** 
*  * 

*  Plots  glint  ratio  versus  receiver  elevation  and  azimuth.  * 

*  * 

*  USES  Sun,  Tilt,  Function,  Rho.  * 

*  Last  revised;  April  26,  1994.  * 

*  * 


************************************************************************ 

c 


Pi 

R2D 

D2R 

Epsilon  = 
Delta  == 
Onem  = 
Onep  = 
Infinity  = 


4 . *ATAN ( 1 . ) 
180. /Pi 
Pi/180. 
D2R*0.2659 
1.4E-6 
1.  -  Delta 
1.  +  Delta 
999999 


Zero  =  0. 

TO  =  273.15 


C 


C 

C 


C 


REAL*4  NoverNo 


IPR  =  10 

OPEN  (IPR,  CARRIAGE  CONTROL  =  'FORTRAN',  FILE  =  'G.DAT') 


PRINT  *,  "Wind  velocity  (m  s-1)  =  " 

READ  (5,  *)  W 
IF  (W  .LE.  0.01)  W  =  0.01 
PRINT  *,  "Solar  Zenith  Angle,  Azimuth  (deg)  =  " 

READ  (5,  *)  To,  Po 
To  =  D2R*To 
Po  =  D2R*Po 

PRINT  *,  "Receiver  Zenith  Angle  (deg):  Min,  Max,  Step  =  " 

READ  (5,  *)  Trmin,  Trmax,  Trstep 
Trmin  =  D2R*Trmin 
Trmax  =  D2R*Trmax 
Trstep  =  D2R*Trstep 

PRINT  *,  "Receiver  Azimuth  (deg);  Min,  Max,  Step  =  " 

READ  (5,  *)  Prmin,  Prmax,  Prstep 
Prmin  =  D2R*Prmin 
Prmax  =  D2R*Prmax 
Prstep  =  D2R*Prstep 

PRINT  *,  "Wavelength  (urn)  for  rho  (enter  0.  for  rho  =  100%)  =  " 
READ  (5,  *)  WL 
IF  (WL  .EQ.  0.)  THEN 

V  =  0. 

ELSE 

V  =  1.E4/WL 

END  IF 

WRITE  (IPR,  10)  "\  A  Sun  Glint  'Picture':  " 
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C; \COX\PLOTGLNT.FOR  4/28/94 


WRITE 

(IPR, 

10) 

"\  N/No  versus 

Pixel  Position" 

WRITE 

(IPR, 

50) 

"\Solar  zenith  (deg)  = 

II 

i 

R2D*To 

WRITE 

(IPR, 

50) 

"\Solar  azimuth  (deg)  = 

II 

/ 

R2D*Po 

WRITE 

(IPR, 

50) 

"\Wavelength  (um)  = 

It 

9 

WL 

WRITE 

(IPR, 

50) 

"\Wind  Speed  (m  s-1)  = 

It 

9 

W 

WRITE 

(IPR, 

150) 

WRITE 

(IPR, 

100) 

ZERO 

DO  Tr  =  Tinnin,  Onep*Tnnax,  Trstep 
WRITE  (IPR,  100)  R2D*Tr 
END  DO 

WRITE  (IPR,  150) 

C 

DO  Pr  =  Prmin,  Onep*Prmax,  Prstep 
WRITE  (IPR,  100)  R2D*Pr 
DO  Tr  =  Trmin,  Onep*Trmax,  Trstep 
CALL  Sun(V,  W,  NoverNo) 

WRITE  (IPR,  200)  NoverNo 
END  DO 

WRITE  (IPR,  150) 

END  DO 
C 

PRINT  *,  "DONE.  DATA  IN  'G.DAT'." 

C 


10  FORMAT 
50  FORMAT 
100  FORMAT 
150  FORMAT 
200  FORMAT 
C 

END 


(A50) 

(/,  A40,  F10.2) 
(F10.2,  \) 

(/) 

(1PE10.2,  \) 
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C:\COX\G.DAT  4/28/94 


\  A  Sun  Glint  'Picture': 

\  N/No  versus  Pixel  Position 

\Solar  zenith  (deg)  =  80.00 

\Solar  azimuth  (deg)  =  90.00 

\Wavelength  (um)  =  0.00 


\Wind 

Speed  (m 

S-1) 

10.00 

0.00 

76.00 

78.00 

80.00 

82.00 

BA  .  00 

260.00 

3.80E-05 

2.82E-05 

1.81E-05 

9.39E-06 

3  •  — ^  '  w  o 

262.00 

8.99E-05 

7.89E-05 

6.29E-05 

4.42E-05 

2.56E-05 

264.00 

1.75E-04 

1.74E-04 

1.65E-04 

1.46E-04 

1.17E-04 

266.00 

2.80E-04 

3.06E-04 

3.26E-04 

3.39E-04 

3.42E-04 

268.00 

3.72E-04 

4.28E-04 

4.91E-04 

5.62E-04 

6.49E-04 

270.00 

4.08E-04 

4.79E-04 

5.62E-04 

6.64E-04 

8.02E-04 

272.00 

3.72E-04 

4.28E-04 

4.91E-04 

5.62E-04 

6.49E-04 

274.00 

2.80E-04 

3.06E-04 

3.26E-04 

3.39E-04 

3.42E-04 

276.00 

1.75E-04 

1.74E-04 

1.65E-04 

1.46E-04 

1.17E-04 

278.00 

8.99E-05 

7.89E-05 

6.29E-05 

4.42E-05 

2.56E-05 

280.00 

3.81E-05 

2.82E-05 

1.81E-05 

9.39E-06 

3.55E-06 
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C:\COX\PLOTP.FOR  4/28/94 


COMMON/Constants/Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

C 

*****************************  PLOTP.FOR  ******************************* 
*  * 

*  Calculates  Cox-Munk  occurence  probability  p.  Prints  p  to  * 

*  a  file  for  subsequent  graphing  on  a  three  dimensional  plot.  * 

*  * 

*  USES  FUNCTION  * 

*  Last  revised:  April  22,  1994  * 

'  *  * 

***********  *****************************11*********  ********  ************* 
C 

=  4.*ATAN(1.) 

=  180. /Pi 
=  Pi/180. 

=  D2R*0.2659 
=  1.4E-6 
=  1.  -  Delta 
=  1.  +  Delta 


Pi 
R2D 
D2R 

Epsilon 
Delta 
Onem 
Onep 

Infinity  =  999999 
Zero  =  0. 

TO  =  273.15 


IPR  =  10 

OPEN  (IPR,  CARRIAGE  CONTROL  =  'FORTRAN’,  FILE  =  'P.DAT') 

PRINT  *,  "CURRENT  WIND  SPEED  (m  s-1)  =  " 

READ  (5,  *)  W 
IF  (W  .LT.  0.01)  THEN 
W  =  0.01 

END  IF 

PRINT  *,  "MINIMUM  X  SLOPE  =  " 

READ  (5,  *)  Sxmin 
PRINT  *,  "MAXIMUM  X  SLOPE  =  " 

READ  (5,  *)  Sxmax 
PRINT  *,  "X  SLOPE  INCREMENT  =  " 

READ  (5,  *)  dSX 
PRINT  *,  "MINIMUM  Y  SLOPE  =  " 

READ  (5,  *)  Symin 
PRINT  *,  "MAXIMUM  Y  SLOPE  =  " 

READ  (5,  *)  Symax 
PRINT  *,  "Y  SLOPE  INCREMENT  =  " 

READ  (5,  *)  dSy 


WRITE  (IPR,  50)  "\OCCURRENCE  PROS.  FOR  WIND  SPEED  =  ", 
+  W,  "M  S-1" 

C 

WRITE  (IPR,  100)  ZERO 
DO  Sy  =  Symin,  Onep*Symax,  dSy 
WRITE  (IPR,  100)  Sy 
END  DO 

WRITE  (IPR,  150) 

C 

DO  Sx  =  Sxmin,  Onep*Sxmax,  dSx 
WRITE  (IPR,  100)  Sx 
DO  Sy  =  Symin,  Onep* Symax,  dSy 
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C;\COX\PLOTP.FOR  4/28/94 


WRITE  (IPR,  200)  p(Sx,  Sy,  W) 
END  DO 

WRITE  (IPR,  150) 

END  DO 
C 

PRINT  *,  "DONE.  DATA  IN  'P.DAT*.'* 

C 

50  FORMAT  (A40,  F6.2,  A6,  //) 

100  FORMAT  (F8.4,\) 

150  FORMAT  (/) 

200  FORMAT  (F8.2,  \) 

C 

END 
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C:\COX\P.DAT  4/28/94 


\OCCURRENCE  PROS.  FOR  WIND  SPEED  -  10.00  MS- 1 


0.0000 

-0.2000 

-0.1000 

0.0000 

0.1000 

0.2000 

0.2000 

1.30 

2.55 

3.19 

2.55 

1.30 

0.1000 

2.08 

4.10 

5.13 

4.10 

2.08 

0.0000 

2.44 

4.80 

6.01 

4.80 

2.44 

0.1000 

2.08 

4.10 

5.13 

4.10 

2.08 

0.2000 

1.30 

2.55 

3.19 

2.55 

1.30 
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C; \COX\PLOTQ.FOR  4/28/94 


COMMON/Constants/Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 


*****************************  PLOTQ. FOR  ******************************* 

*  * 

*  Calculates  Cox-Munk  interaction  probability  q.  Prints  q  to  * 

*  a  file  for  subsequent  graphing  on  a  three  dimensional  plot.  * 

*  * 

*  USES  FUNCTION  * 

*  Last  revised:  April  4,  1994  * 

K*  * 

*********************************************************************** 


Pi  =  4.*ATAN(1.) 

R2D  =  180. /Pi 

D2R  =  Pi/180. 

Epsilon  =  D2R*0.2659 
Delta  =  1.4E-6 
Onem  =  1.  -  Delta 
Onep  =  1.  +  Delta 
Infinity  =  999999 
Zero  =  0. 

TO  =  273.15 
N  =2 
M  =7 

IPR  =  10 

OPEN  (IPR,  CARRIAGE  CONTROL  =  'FORTRAN',  FILE  =  'Q.DAT') 

PRINT  *,  "CURRENT  WIND  SPEED  (m  s-1)  =  " 

READ  (5,  *)  W 
IF  (W  .LT.  0.01)  THEN 
W  =  0.01 

END  IF 

PRINT  *,  "BEAM  ZENITH  ANGLE  (deg)  =  " 

READ  (5,  *)  T 
T  =  D2R*T 

PRINT  *,  "BEAM  AZIMUTH  ANGLE  (deg)  =  " 

READ  (5,  *)  P 
P  =  D2R*P 

PRINT  *,  "MINIMUM  X  SLOPE  =  " 

READ  (5,  *)  Sxmin 
PRINT  *,  "MAXIMUM  X  SLOPE  =  " 

READ  (5,  *)  Sxmax 
PRINT  *,  "X  SLOPE  INCREMENT  =  " 

READ  (5,  *)  dSx 
PRINT  *,  "MINIMUM  Y  SLOPE  =  " 

READ  (5,  *)  Symin 
PRINT  *,  "MAXIMUM  Y  SLOPE  =  " 

READ  (5,  *)  Symax 
PRINT  *,  "Y  SLOPE  INCREMENT  =  " 

READ  (5,  *)  dSy 

WRITE  (IPR,  50)  "\ INTERACTION  PROB.  FOR  WIND  SPEED  =  ", 

+  W,  "M  S-1" 

WRITE  (IPR,  50)  "\ ZENITH  ANGLE  OF  BEAM  =  ", 

+  R2D*T,  "DEG" 
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WRITE  (IPR,  50)  "\AZIMUTHAL  ANGLE  OF  BEAM  =  , 

+  R2D*P,  "DEG" 

WRITE  (IPR,  60)  "N INTEGRAL  =  ,  BoverA(T, P, W) 

WRITE  (IPR,  60)  'VJOS  (ZENITH  ANGLE)  =  COS(T) 

WRITE  (IPR,  60)  *'\ INTEGRAL/ COS  =  ”,  BoverA(T,P,W) /COS(T) 
WRITE  (IPR,  70)  ”\N  =  ”,  N 
WRITE  (IPR,  70)  ”\M  =  ”,  M 

WRITE  (IPR,  100)  ZERO 
DO  Sy  =  Symin,  Symax,  dSy 
WRITE  (IPR,  100)  Sy 
END  DO 

WRITE  (IPR,  150) 

DO  Sx  =  Sxmln,  Sxmax,  dSx 
WRITE  (IPR,  100)  SX 
DO  Sy  =  symin,  Symax,  dSy 

WRITE  (IPR,  200)  q(T,  P,  Sx,  Sy,  W) 

END  DO 

WRITE  (IPR,  150) 

END  DO 

PRINT  *,  "DONE.  DATA  IN  'Q.DAT*.” 

50  FORMAT  (A40,  F6,2,  A6,  //) 

60  FORMAT  (A40,  F9.5,  //) 

70  FORMAT  (A40,  13,  //) 

100  FORMAT  (F8.4,\) 

150  FORMAT  (/) 

200  FORMAT  (F8.2,  \) 

END 
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C; \COX\Q.DAT  4/28/94 


\INTERACTION  PROB.  FOR  WIND  SPEED  =  10.00  M  S-1 

\ZENITH  ANGLE  OF  BEAM  «  80.00  DEG 

\AZIMUTHAL  ANGLE  OF  BEAM  =  270.00  DEG 

\INTEGRAL  =  0.18235 

\COS  (ZENITH  ANGLE)  *  0.17365 

\INTEGRAL/COS  *  1.05014 


\N  =  2 
\M  =  7 


0.0000 

-0.2000 

-0.1000 

0.0000 

0.1000 

0.2000 

0.2000 

0.00 

1.05 

3.04 

3.80 

2.63 

0.1000 

0.00 

1.69 

4.88 

6.11 

4.23 

0.0000 

0.00 

1.98 

5.72 

7.16 

4.96 

0.1000 

0.00 

1.69 

4.88 

6.11 

4.23 

0.2000 

0.00 

1.05 

3.04 

3.80 

2.63 
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C:\COX\PLOTRCVR.FOR  4/28/94 


C 

COMMON/ IN/TO , Po , Tr , Pr 
COMMON/OUT/Sx, Sy , Omega , Tn , Pn 
COMMON/Constants/Pl,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

C 

*****************************  piotRcvr . FOR  ***************************** 
*  * 

*  To  plot  glint  as  a  function  of  zenith  angle  of  receiver.  * 

*  * 

'*  USES  Sun,  Tilt,  Function,  Rho.  * 

*  Latest  revision:  April  26,  1994  * 

*  * 

************************************************************************ 

’c 

Pi  =  4 . *ATAN ( 1 . ) 

R2D  =  180. /Pi 

D2R  =  Pi/180. 

Epsilon  =  D2R*0.2659 
Delta  =  1.4E-6 
Onem  =  1.  -  Delta 
Onep  =  1.  +  Delta 
Infinity  =  999999 
Zero  =  0. 

TO  =  273.15 
C 

REAL*4  NoverNo 
C 

IPR  =  10 

OPEN  (IPR,  CARRIAGE  CONTROL  =  'FORTRAN*,  FILE  =  'R.DAT') 

C 

PRINT  *,  "Solar  Zenith  Angle,  Azimuth  (deg)  =  " 

READ  (5,  *)  To,  Po 
To  =  D2R*To 
Po  =  D2R*Po 

PRINT  *,  "Receiver  Zenith  Angle  (deg):  Min,  Max,  Step  =  " 

READ  (5,  *)  Trmin,  Trmax,  Trstep 
Trmin  =  D2R*Trmin 
Trmax  =  D2R*Trmax 
Trstep  =  D2R*Trstep 
PRINT  *,  "Receiver  Azimuth  (deg)  =  " 

READ  (5,  *)  Pr 
Pr  =  D2R*Pr 

PRINT  *,  "Wind  Speed  (m/s)  =  " 

READ  (5,  *)  W 

IF  (W  .EQ.  0.)  W  =  0.01 

PRINT  *,  "Wavelength  (urn)  for  rho  (enter  0.  for  rho  =  100%)  =  " 
READ  (5,  *)  WL 
IF  (WL  .EQ.  0.)  THEN 

V  =  0. 

ELSE 

V  =  1.E4/WL 

END  IF 
C 

WRITE  (IPR,  50)  "\Solar  zen,  az  (deg)  =  ",  R2D*To,  R2D*Po 

WRITE  (IPR,  60)  "\Receiver  az  (deg)  =  ",  R2D*Pr 

WRITE  (IPR,  50)  "\W  (m  s-1) ,  Lambda  (urn)  =  ",  W,  WL 
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WRITE  (IPR,  75) 

+''\  Tr  (deg)  Nsun/No  Flat-top  Cox-Munk  " 

C 

DO  Tr  =  Trillin,  Trmax,  Trstep 
CALL  TILT 

C  at  the  solar  center  to  calculate  the  Cox-Munk  glint  ratio: 

Gamma  *  l./(4.*COS(Tr)*COS(Tn)**4) 

G  =  Gamma*Pi*Epsilon**2 

RpG  ==  Rho(Omega,  V)*p(Sx,  Sy,  W)  *G 

C  and  to  calculate  the  flat-top  glint  ratio; 

Sigma  =  1. / (4 . *COS(Omega) *COS(Tn) **3) 

S  =»  Sigma*Pi*Epsilon**2 

RqS  =  Rho(Omega,  V)*q(Tr,  Pr,  Sx,  Sy,  W) *S 
C  and  then  get  the  actual  glint  ratio: 

CALL  Sun(V,  W,  NoverNo) 

WRITE  (IPR,  100)  R2D*Tr,  NoverNo,  RqS,  RpG 
END  DO 
C 

50  FORMAT  (A35,2F10.2) 

60  FORMAT  (A35,F20.2) 

75  FORMAT  (//,  A58) 

100  FORMAT  (F12.2,  3(1PE15.4)) 

C 

PRINT  *,  "DONE.  DATA  IN  'R.DAT*." 

C 

END 


(39) 

(24) 

(25) 
(35) 
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C:\COX\R.DAT  4/28/94 


\Solar  zen, 

az  (deg)  == 

80.00 

90.00 

\Receiver  az 

i  (deg) 

270.00 

\W  (m  s-l) , 

Lambda  (um)  = 

10.00 

0.00 

Tr  (deg) 

Nsun/No 

Flat-top 

Cox-Munk 

70.00 

2.5603E-04 

2.5397E-04 

2.5397E-04 

72.00 

2.9992E-04 

2.9751E-04 

2.9751E-04 

74.00 

3.4941E-04 

3.4663E-04 

3.4855E-04 

76.00 

4.0810E-04 

4.0488E-04 

4.0977E-04 

78.00 

4.7906E-04 

4.7532E-04 

4.8583E-04 

80.00 

5.6172E-04 

5.5740E-04 

5.8534E-04 

82.00 

6.6431E-04 

6.5930E-04 

7.2579E-04 

84.00 

8.0166E-04 

7.9577E-04 

9.4837E-04 

86.00 

9.6035E--04 

9.5358E-04 

1.3773E-03 

88.00 

1.1692E-03 

1.1614E-03 

2.6344E-03 

90.00 

1.4967E-03 

1.4879E-03 

-2.1643E+02 
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C;\COX\PLOTSUN.FOR  4/28/94 


C 

COMMON/ IN/To , Po , Tr , Pr 
COMMON/OUT/  Sx ,  Sy ,  Onega ,  Tn ,  Pn 
COMMON/Constants/Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onen,  Onep,  Infinity,  Zero,  TO 

C 

*****************************  piotSun. FOR  ****************************** 


*  * 

*  To  plot  glint  as  a  function  of  solar  zenith  angle.  * 

*  * 

*  USES  Sun,  Tilt,  Function,  Rho.  * 

*  Latest  revision:  April  26,  1994  * 

*  * 


************************************************************************ 

c 


Pi 

R2D 

D2R 

Epsilon 

Delta 

Onen 

Onep 

Infinity 


4.*ATAN(1.) 
180. /Pi 
Pi/180. 
D2R*0.2659 
1.4E-6 
1 .  -  Delta 
1 .  +  Delta 
999999 


Zero  =  0. 

TO  =  273.15 


C 


C 


REAL*4  NoverNo 


C 


C 


IPR  =  10 

OPEN  (IPR,  CARRIAGE  CONTROL  =  'FORTRAN’,  FILE  =  'S.DAT') 


PRINT  *,  "Solar  Zenith  Angle  (deg):  Min,  Max,  Step  =  " 

READ  (5,  *)  Tonin,  Tonax,  Tostep 
Tonin  =  D2R*Tonin 
Tonax  =  D2R*Tonax 
Tostep  =  D2R*Tostep 

PRINT  *,  "Solar  Azinuth  (deg)  =  " 

READ  (5,  *)  Po 
Po  =  D2R*Po 

PRINT  *,  "Receiver  Zenith  Angle,  Azinuth  (deg)  =  " 

READ  (5,  *)  Tr,  Pr 
Tr  =  D2R*Tr 
Pr  =  D2R*Pr 

PRINT  *,  "Wind  Speed  (n/s)  =  " 

READ  (5,  *)  W 

IF  (W  .EQ.  0.)  W  =  0.01 

PRINT  *,  "Wavelength  (un)  for  rho  (enter  0.  for  rho  =  100%)  =  " 
READ  (5,  *)  WL 
IF  (WL  .EQ.  0.)  THEN 

V  =  0. 

ELSE 

V  =  1.E4/WL 

END  IF 


WRITE  (IPR,  50)  "\Receiver  zen,  az  (deg) 
WRITE  (IPR,  60)  "\Solar  az  (deg) 

WRITE  (IPR,  50)  "\W  (n  s-1).  Lambda  (un) 


",  R2D*Tr,  R2D*Pr 
",  R2D*Po 
",  W,  WL 
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WRITE  (IPR,  75) 

+''\  To  (deg)  Nsun/No  Flat-top  Cox-Munk  ” 

C 

DO  To  =  Tomin,  Tomax,  Tostep 
CALL  TILT 

C  at  the  solar  center  to  calculate  the  Cox-Munk  glint  ratio: 


Gamma  =  1. / (4 . *COS(Tr) *COS(Tn) **4) 

G  =  Gamma*Pi*Epsilon**2  (39) 

RpG  =  Rho(Omega,  V)*p(Sx,  Sy,  W)*G 

C  and  to  calculate  the  flat-top  glint  ratio: 

Sigma  =  1. / (4 . *COS (Omega) *COS(Tn) **3)  (24) 

S  =  Sigma*Pi*Epsilon**2  (25) 

RqS  =  Rho(Omega,  V)*q(Tr,  Pr,  Sx,  Sy,  W) *S  (35) 


C  and  then  get  the  actual  glint  ratio: 

CALL  Sun(V,  W,  NoverNo) 

WRITE  (IPR,  100)  R2D*To,  NoverNo,  RqS,  RpG 
END  DO 
C 

50  FORMAT  (A35,2F10.2) 

60  FORMAT  (A35,F20.2) 

75  FORMAT  (//,  A58) 

100  FORMAT  (F12.2,  3(1PE15.4)) 

C 

PRINT  *,  "DONE.  DATA  IN  'S.DAT'." 

C 

END 
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C; \COX\S.DAT  4/28/94 


\Receiver  zen,  az 

(deg)  =  80. 

00 

270.00 

\Solar  az  (deg) 

= 

90.00 

\W  (m  s-1) ,  Lambda 

(um)  =  10. 

00 

0.00 

To  (deg) 

Nsun/No  Flat-top 

Cox-Munk 

70.00 

4.8019E-04 

4.7634E-04 

5.0022E-04 

71.00 

4.9480E-04 

4.9084E-04 

5.1545E-04 

72.00 

5.0822E-04 

5.0417E-04 

5.2944E-04 

73.00 

5.2034E-04 

5.1620E-04 

5.4208E-04 

74.00 

5.3106E-04 

5.2684E-04 

5.5326E-04 

75.00 

5.4028E-04 

5.3601E-04 

5.6288E-04 

76.00 

5.4793E-04 

5.4362E-04 

5.7088E-04 

77.00 

5.5395E-04 

5.4961E-04 

5.7716E-04 

78.00 

5.5828E-04 

5.5392E-04 

5.8170E-04 

79.00 

5.6087E-04 

5.5653E-04 

5.8443E-04 

80.00 

5.6172E-04 

5.5740E-04 

5.8534E-04 

81.00 

5.6081E-04 

5.5653E-04 

5.8443E-04 

82.00 

5.5815E-04 

5.5392E-04 

5.8170E-04 

83.00 

5.5375E-04 

5.49eiE-04 

5.7716E-04 

84.00 

5.4766E-04 

5.43o2E-04 

5.7088E-04 

85.00 

5.3993E-04 

5.3601E-04 

5.6288E-04 

86.00 

5.3061E-04 

5.2684E-04 

5.5326E-04 

87.00 

5.1979E-04 

5.1620E-04 

5.4208E-04 

88.00 

5.0755E-04 

5.0417E-04 

5.2944E-04 

89.00 

4.9399E-04 

4.9084E-04 

5.1545E-04 

90.00 

4.7886E-04 

4.7634E-04 

5.0022E-04 
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C 

COMMON/ IN/To , Po , Tr , Pr 

COMMON/OUT/ Sx , Sy , Omega , Tn , Pn 

COMMON/ Constants /Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

C 

it  it  it  it  "k  it  it  it  it  it  it  ii'k'k’k  it  it  it  it  it  it  it  it  it  it  is  it  it  it  PlO'tWi.nd  •  FOR  ititititititititititititititititifitititititititititititit 
it  * 

*  To  plot  glint  as  a  function  of  wind  velocity.  * 

*  * 

*  USES  Sun,  Tilt,  Function,  Rho.  * 

*  Latest  revision:  April  26,  1994  * 

*  * 

c 

Pi  =  4 . *ATAN ( 1 . ) 

R2D  =  180. /Pi 

D2R  =  Pi/180. 

Epsilon  =  D2R*0.2659 
Delta  =  1.4E-6 
Onem  =  1.  -  Delta 
Onep  =  1.  +  Delta 
Infinity  =  999999 
Zero  =  0. 

TO  =  273.15 
C 

REAL* 4  NoverNo 
C 

IPR  *=  10 

OPEN  (IPR,  CARRIAGE  CONTROL  =  'FORTRTU^',  FILE  =  'W.DAT') 

C 

PRINT  *,  "Solar  Zenith  Angle,  Azimuth  (deg)  =  " 

READ  (5,  *)  TO,  Po 
To  =  D2R*To 
Po  =  D2R*Po 

PRINT  *,  "Receiver  Zenith  Angle,  Azimuth  (deg)  =  " 

READ  (5,  *)  Tr,  Pr 
Tr  =  D2R*Tr 
Pr  =  D2R*Pr 

PRINT  *,  "Wavelength  (urn)  for  rho  (enter  0.  for  rho  =  100%)  =  " 
READ  (5,  *)  WL 
IF  (WL  .EQ.  0.)  THEN 

V  =  0. 

ELSE 

V  =  1.E4/WL 

END  IF 

PRINT  *,  "Wind  velocity  (m  s-1) ;  Min,  Max,  Step  =  " 

READ  (5,  *)  Wmin,  Wmax,  Wstep 
IF  (Wmin  .EQ.  Zero)  Wmin  =  0.01 
C 

WRITE  (IPR,  50)  "\Solar  zen,  az  (deg)  =  ",  R2D*To,  R2D*Po 

WRITE  (IPR,  50)  "\Recvr  zen,  az  (deg)  =  ",  R2D*Tr,  R2D*Pr 

WRITE  (IPR,  60)  "\Lambda  (urn)  =  ",  WL 

WRITE  (IPR,  75) 

+"\  W  (m/s)  Ng/No  Flat  Top  Cox-Munk  " 

C 

DO  W  =  Wmin,  Wmax,  Wstep 
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CALL  TILT 

C  at  the  solar  center  to  calculate  the  Cox-Munk  glint  ratio: 

Gamma  =  1./ (4.*COS(Tr) *COS(Tn) **4) 

G  =  Gamma*Pi*Epsilon**2 
RpG  =  Rho(Omega,  V)*p(Sx,  Sy,  W)*G 
C  and  to  calculate  the  approximate  glint  ratio: 

Sigma  =  l. / (4 . *COS(Omega) *COS(Tn) **3) 

S  =  Sigma*Pi*Epsilon**2 

RqS  =  Rho (Omega,  V)*q(Tr,  Pr,  Sx,  Sy,  W) *S 
C  and  then  get  the  actual  glint  ratio: 

CALL  Sun(V,  W,  NoverNo) 

WRITE  (IPR,  100)  W,  NoverNo,  RqS,  RpG 
END  DO 
C 

50  FORMAT  (A35,2F10.2) 

60  FORMAT  (A35,F10.2) 

75  FORMAT  (//,  A58) 

100  FORMAT  (F12.2,  3(1PE15.4)) 

C 

PRINT  * ,  '•  DONE .  DATA  IN  '  W .  DAT  •  . " 

C 

END 


(39) 

(24) 

(25) 
(35) 
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C:\COX\W.DAT  4/28/94 


\Solar  zen,  az 

(deg)  = 

80.00 

90.00 

\Recvr  zen,  az 

(deg)  = 

80.00 

270.00 

\ Lambda  (um) 

= 

0.00 

W  (m/s) 

Ng/No 

Flat  Top 

Cox-Munk 

1.00 

3.9374E-03 

3.9319E-03 

3.9319E-03 

2.00 

2.3696E-03 

2.3580E-03 

2.3580E-03 

3.00 

1.7015E-03 

1.6913E-03 

1.7013E-03 

4.00 

1.3283E-03 

1.3195E-03 

1.3344E-03 

5.00 

1.0889E-03 

1.0813E-03 

1.0988E-03 

6.00 

9.2206E-04 

9.1540E-04 

9.3439E-04 

7.00 

7.9852E-04 

7.9262E-04 

8.1299E-04 

8.00 

7.0138E-04 

6.9611E-04 

7.1962E-04 

9.00 

6.2460E-04 

6.1984E-04 

6.4555E-04 

10.00 

5.6172E-04 

5.5740E-04 

5.8534E-04 

11.00 

5.0975E-04 

5.0580E-04 

5.3543E-04 

12.00 

4.6618E-04 

4.6254E-04 

4.9337E-04 

13.00 

4.2916E-04 

4.2578E-04 

4.5745E-04 

14.00 

3.9731E-04 

3.9417E-04 

4.2642E-04 

15.00 

3.6963E-04 

3.6670E-04 

3.9933E-04 

16.00 

3.4532E-04 

3.4257E-04 

3.7548E-04 

17.00 

3.2373E-04 

3.2115E-04 

3.5432E-04 

18.00 

3.0453E-04 

3.0209E-04 

3.3542E-04 

19.00 

2.8735E-04 

2.8504E-04 

3.1844E-04 

20.00 

2.7189E-04 

2.6970E-04 

3.0309E-04 
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C:\COX\TESTB.FOR  1I1AI9A 


COMMON/ Constants /Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

C 

TestB.  FOR  Itififklfklfklfkifkltiftfk'klfkltlflfk'ifk-k-k-k'k'k 
*  * 

*  To  display  a  value  of  BoverA  by  entering  arguments  at  the  * 

*  keyboard .  * 

•k  k 

*  USES  Function.  * 

*  Last  revised  June  6,  1994.  * 


C 

Pi  =  4.*ATAN(1.) 

R2D  =  180. /Pi 

D2R  =  Pi/180. 

Epsilon  =  D2R*0.2659 
Delta  =  1.4E-6 
Onem  =  1.  -  Delta 
Onep  =  1.  +  Delta 
Infinity  =  999999 
Zero  =  0. 

TO  =  273.15 

c 

PRINT  *,  "Beam  zenith  angle,  azimuth  (deg)  =  " 
READ  (5,  *)  T,  P 
T  =  D2R*T 
P  =  D2R*P 

PRINT  *,  "Wind  speed  =  " 

READ  (5,  *)  W 
IF  (W  .LT.  0.01)  W  =  0.01 
C 


C 


PRINT  *,  "BoverA  =  ",  BoverA (T,  P,  W) 


END 
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C:\COX\TESTP.FOR  4/28/94 


COMMON/ Constants /Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

C 

*******************************  TESTP. FOR  ****************************** 
*  * 

*  To  display  a  value  of  p  by  entering  arguments  at  the  keyboard.  * 

*  * 

*  USES  Function.  * 

*  Last  revised  16  March,  1994.  * 

*  * 
************************************************************************ 
C 

Pi  =  4.*ATAN(1.) 

R2D  =  180. /Pi 

D2R  =  Pi/180. 

Epsilon  =  D2R*0.2659 

Delta  =  1.4E-6 

Onem  =  1.  -  Delta 

Onep  =  1.  +  Delta 

Infinity  =  999999 

Zero  =  0 . 

TO  =  273.15 

C 

PRINT  *,  "Sx,  Sy  =  •• 

READ  (5,  *)  Sx,  Sy 
PRINT  *,  "Wind  speed  =  '• 

READ  (5,  *)  W 
IF  (W  .LT.  0.01)  W  =  0.01 
C 

PRINT  *,  "p  =  p(Sx,  Sy,  W) 

C 

END 


Page  1 


C; \COX\TESTQ.FOR  4/28/94 


COMMON/ Constants /Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

C 

*******************************  TestQ.FOR  ****************************** 

*  * 

*  To  display  a  value  of  q  by  entering  argviments  at  the  keyboard.  * 

*  * 

*  USES  Function  * 

*  Last  revised:  April  4,  1994  * 

*  * 
************************************************************************ 


C 


C 

C 


PRINT  *,  "Source  zenith,  azimuth  =  " 
READ  (5,  *)  Ts,  Ps 
Ts  =  D2R*Ts 
Ps  =  D2R*Ps 
PRINT  *,  "Sx,  Sy  =  " 

READ  (5,  *)  Sx,  Sy 
PRINT  *,  "Wind  speed  =  " 

READ  (5,  *)  W 

IF  (W  .LT.  0.01)  W  =  0.01 

PRINT  *,  "q  =  ",  q(Ts,  Ps,  Sx,  Sy,  W) 

END 
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C:\COX\TESTILT.FOR  4/28/94 


C 

COMMON/ IN/Ths , Phs , Thr , Phr 

COMMON/OUT/ Sx , Sy , Omega , Tn , Pn 

COMMON/ Constants /Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

C 

*****************************  TesTilt .FOR  ****************************** 
*  * 

*  To  test  Tilt. FOR  at  the  keyboard.  * 

*  * 

*  USES  Tilt  * 

*  Latest  revision:  March  11,  1994  * 

*  * 
************************************************************************ 
C 

Pi  =  4.*ATAN(1.) 

R2D  =  180. /Pi 

D2R  =  Pi/180. 

Epsilon  =  D2R*0.2659 
Delta  =  1.4E-6 
Onem  =  1.  -  Delta 
Onep  =  1.  +  Delta 
Infinity  =  999999 
Zero  =  0. 

TO  =  273.15 
C 

PRINT  *,  "Source  Zenith  Angle,  Azimuth  (deg)  " 

READ  (5,  *)  Ts,  PS 
Ths  =  D2R*Ts 
Phs  =  D2R*Ps 
C 

PRINT  *,  "Receiver  Zenith  Angle,  Azimuth  (deg)  " 

READ  (5,  *)  Tr,  Pr 
Thr  =  D2R*Tr 
Phr  =  D2R*Pr 
C 

CALL  Tilt 
C 

PRINT  *,  "Sx  =  ",  Sx,  "(Alpha  =  ",  R2D*ATAN(SX) ,  "  deg)" 

PRINT  *,  "-tan*cos  =  ",  -TAN(Tn) *COS(Pn) 

PRINT  * 

PRINT  *,  "Sy  =  ",  Sy,  "(Beta  =  ",  R2D*ATAN(Sy) ,  "  deg)" 

PRINT  *,  "-tan*sin  =  ",  -TAN(Tn) *SIN(Pn) 

PRINT  * 

PRINT  *,  "Omega  =  ",  R2D*Omega,  "  deg" 

PRINT  * 

PRINT  *,  "Tilt  =  ",  R2D*Tn,  "  deg" 

PRINT  *,  "Tilt  Az.  =  ",  R2D*Pn,"  deg" 

C 

END 
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C:\COX\TESTSEA.FOR  IIIAISA 


C 

COMMON/ IN/  Ts,Ps,Tr,Pr 

COMMON/OUT/  Sx,Sy, Omega, Tn,Pn 

COMMON /Constants/  Pi , R20, D2R, Epsilon, 

+  Delta , Onem, Onep, Infinity, Zero, TO 

COMMON/ Sea/  a,b,c,Tsea 
C 

*****************************  TestSea. FOR  ****************************** 
*  * 


*  To  test  Sky. FOR  and  Sun. FOR  from  the  keyboard.  * 

it  it 

*  USES  Sun,  Sky,  Tilt,  Function,  Rho,  and  an  input  file  Sea. In.  * 

*  Latest  revision:  April  26,  1994  * 

*  * 


************************************************************************ 


C 


C 

C 


C 


C 


C 


C 


Pi 

= 

4 . *ATAN ( 1 . ) 

R2D 

180. /Pi 

D2R 

=s 

Pi/180. 

Epsilon 

s= 

D2R*0.2659 

Delta 

= 

1.4E-6 

Onem 

= 

1 .  -  Delta 

Onep 

s 

1 .  +  Delta 

Infinity 

SR 

999999 

Zero 

= 

0. 

TO 

273.15 

Tsun 

= 

5900. 

REAL* 4  Ns 

’  f 

No,  Nsky,  Nsea 

IRD  =  15 

OPEN  (IRD, 

FILE  =  'SEA. IN 

PRINT  *, 

"Receiver  Zenith 

READ 

(5 

i,  *)  Tr,  Pr 

Tr  =  D2R*Tr 
Pr  =  D2R*Pr 


PRINT  *,  "Wind  Speed  (m/s)  =  " 
READ  (5,  *)  W 
IF  (W  .EQ.  0.)  W  =  0.01 


•• 


READ  (IRD,  100)  V 

READ  (IRD,  200)  Tso,  Pso,  Tau 

Ts  =  D2R*Tso 
Ps  =  D2R*Pso 
READ  (IRD,  300)  a,  b,  C 

READ  (IRD,  100)  Tsea 

Tsea  =  Tsea  +  TO 


CALL  Sun(V,  W,  NoverNo) 

No  =  BB(V,  Tsun) *Tau 

Nsun  =  NoverNo*No 

CALL  Sky(V,  W,  Nsky,  Nsea) 

PRINT  *,  "Under  the  following  conditions:" 

PRINT  *,  "  Wind  Speed  (m/s)  =  ",  W 

PRINT  *,  "  Receiver  zen.  (deg)  =  ",  R2D*Tr 
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C;\COX\TESTSEA.FOR  7/14/94 


PRINT 

*  ^ 

II 

Receiver  azm.  (deg) 

=  ",  R2D*Pr 

PRINT 

* 

II 

Solar  zenith  (deg) 

=  ",  Tso 

PRINT 

■k 

II 

Solar  azimuth  (deg) 

=  ”,  Pso 

PRINT 

■k 

II 

Wave  Number  ((cm-1) 

=  ",  V 

PRINT 

k  ^ 

II 

Atmos.  Transmission 

=  ",  Tau 

PRINT 

k 

PRINT 

k 

II 

No  at  top  atm. 

BB(V,  Tsun) 

PRINT 

k 

It 

No  at  footprint 

No 

PRINT 

k 

II 

Ns  at  zenith 

=  '• 

Ns(a,b,c,0. 

PRINT 

k  ^ 

II 

Ns  at  80.  deg 

—  "  / 

Ns(a,b,c,80 

PRINT 

k  ^ 

II 

Ns  at  horizon 

=  , 

Ns(a,b,c,90 

PRINT 

k 

II 

N*  for  sea 

BB(V,Tsea) , 

+ 

••  at  T(C) 

Tsea-TO 

PRINT 

k 

PRINT 

k 

"the  Cox-Munk  radiance 

values  are;” 

PRINT 

k  ^ 

II 

Nsun 

_  II 
“  # 

Nsun 

PRINT 

k  ^ 

II 

Nsky 

^  II 

“  9 

Nsky 

PRINT 

k 

II 

Nsea 

—  II 

Nsea 

PRINT 

k 

PRINT 

"Note: 

All  radiance  values 

in  (W  m-2  sr 

100 

FORMAT 

(F10.3) 

200 

FORMAT  (3F10.3) 

300 

FORMAT  (3E10.3) 

END 
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C:\C0X\SEA1.IN 


945.  V  (cm-1)  FORMAT 

80.  90.  0.3939  Ts  (deg),  Ps  (deg),  Tau  FORMAT 

1.1779E-2  5.8410E-5  2.2542E-1  a,  b,  c  (W  in-2  sr-1  (cm-l)-l)  FORMAT 

15.  Tsea  (C)  FORMAT 

This  is  file  "Seal. In"  (Copy  to  "Sea. In") 

With  Receiver  zenith  angle  =  80  deg. 

Receiver  azimuth  =  270  deg. 

Wind  speed  =  10  m  s-l 

the  output  of  "TestSea"  should  be 

Nsun  =  2.51138E-3  W  m-2  sr-1  (cm-l)-l 
Nsky  =  6.31824E-3  W  m-2  sr-1  (cm-l)-l 
Nsea  =  7.09113E-2  W  m-2  sr-1  (cm-l)-l 


F10.3 

3F10.3 

3E10.3 

F10.3 
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C:\COX\TESTSKV.FOR  1I1AI9A 


C 

COMMON/ IN/  Ts,Ps,Tr,Pr 

COMMON/OUT/  Sx , Sy , Omega , Tn , Pn 

COMMON /Constants/  Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

COMMON/Sea/  a,  b,  c,  Tsea 
C 

*****************************  TestSky . FOR  ****************************** 
*  * 

*  To  test  Sky. FOR  from  the  keyboard.  * 

*  * 

*  USES  Sky,  Tilt,  Function,  Rho.  * 

*  Latest  revision:  May  2,  1994  * 

it  it 


it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  ±  it  it  it  it  it  it  it  it  it  it  it  it  it 


c 


c 

c 


c 


c 

c 


Pi 

4.*ATAN(1.) 

R2D 

= 

180. /Pi 

D2R 

= 

Pi/180. 

Epsilon 

= 

D2R*0.2659 

Delta 

= 

1.4E-6 

Onem 

1 .  -  Delta 

Onep 

= 

1.  +  Delta 

Infinity 

= 

999999 

Zero 

= 

0. 

TO 

273.15 

REAL*4  Ns, 

Nsky,  Nsea 

a  =  1.1779E-2 
b  =  5.8410E-5 
C  =  12.915652*D2R 
Tsea  =  TO  +  15. 


PRINT  *,  "Receiver  Zenith  Angle,  Azimuth  (deg)  =  ” 

READ  (5,  *)  Tr,  Pr 
Tr  =  D2R*Tr 
Pr  =  D2R*Pr 

PRINT  *,  "Wind  Speed  (m/s)  =  " 

READ  (5,  *)  W 

IF  (W  .EQ.  0.)  W  =  0.01 

PRINT  *,  "Wavelength  (urn)  for  rho  (enter  0.  for  rho  -  100%)  =  " 
READ  (5,  *)  WL 
IF  (WL  .EQ.  0.)  THEN 

V  =  0. 

ELSE 

V  =  1.E4/WL 

END  IF 


CALL  Sky ( V , W , Nsky . Nsea ) 


PRINT 

PRINT 

* 

"Ns  at  zenith 

=  tl 

t 

Ns(a,b,c,0.) 

PRINT 

"Ns  at  horizon 

~~  11 

t 

Ns(a,b,c,90.*D2R) 

PRINT 

"N*  for  sea 

SS  If 

9 

BB(V,Tsea),  "  at  T(K)  =  ",  Tsea 

PRINT 

PRINT 

* 

"Nsky 

ss  tt 

9 

Nsky 
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C:\COX\TESTSKY.FOR  7/14/94 
PRINT  *,  "Nsea  =  ,  Nsea 

END 
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C:\COX\TESTSUN.FOR  4/28/94 


C 

COMMON/ IN/Ts , Ps , Tr , Pr 

COMMON / OUT /Sx , Sy , Omega , Tn , Pn 

COMMON/ Constants /Pi,  R2D,  D2R,  Epsilon, 

+  Delta,  Onem,  Onep,  Infinity,  Zero,  TO 

C 

Ithlelclfklfk'klfk'klili-kltlelfk'itle'kltlilfklfklt  TeStSun .  FOR  ‘k'kifkltlt'itlfkitlfkitlfifkltlilfkifklilcifk'klilfk 


*  * 

*  To  test  Sun. FOR  from  the  keyboard.  * 

*  * 

*  USES  Sun,  Tilt,  Function,  Rho.  * 

*  Latest  revision:  April  26,  1994  * 

*  * 


iticicicicicicitie'k'kieitic'k'kic'kieic’k'k’kit'k’k'k’kicitic'kitititic’kit’k'kic’k'kieit'k'kit'kic'kieieic'k'k'kicicicieiticicicic'kic’kiciit'it 

'c 

Pi 
R2D 
D2R 

Epsilon 
Delta 
Onem 
Onep 

Infinity 
Zero 
TO 
C 

REAL* 4  NoverNo 

c 

PRINT  *,  "Solar  Zenith  Angle,  Azimuth  (deg)  =  " 

READ  (5,  *)  Ts,  Ps 
Ts  «  D2R*Ts 
Ps  =  D2R*Ps 

PRINT  *,  "Receiver  Zenith  Angle,  Azimuth  (deg)  =  " 

READ  (5,  *)  Tr,  Pr 
Tr  =  D2R*Tr 
Pr  =  D2R*Pr 

PRINT  *,  "Wind  Speed  (m/s)  =  " 

READ  (5,  *)  W 
IF  (W  .EQ.  0.)  W  =  0.01 

PRINT  *,  "Wavelength  (urn)  for  rho  (enter  0.  for  rho  =  100%)  =  " 
READ  (5,  *)  WL 
IF  (WL  .EQ.  0.)  THEN 

V  =  0. 

ELSE 

V  =  1.E4/WL 

END  IF 
C 

CALL  TILT 

C  at  the  solar  center  to  calculate  the  Cox-Munk  glint  ratio: 

Gamma  =  1. / (4 . *COS(Tr) *COS (Tn) **4) 

G  =  Gamma*Pi*Epsilon**2 

RpG  =  Rho(Omega,  V)*p(Sx,Sy,  W)*G 

C  then  to  calculate  the  "flat-top"  glint  ratio: 

Sigma  =  1. / (4 . *COS (Omega) *COS(Tn) **3) 

S  =  Sigma*Pi*Epsilon**2 

RqS  =  Rho(Omega,  V)*q(Tr,  Pr,  Sx,  Sy,  W)*S 

C  and  finally  to  get  the  exact  glint  ratio: 


(39) 

(39) 

(24) 

(25) 
(35) 


=  4.*ATAN(1.) 
=  180. /Pi 
=  Pi/180. 

=  D2R*0.2659 
=  1.4E-6 
=  1.  -  Delta 
=  1.  +  Delta 
=  999999 
=  0. 

=  273.15 
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C:\COX\TESTSUN.FOR  4/28/94 


CALL  Sun(V,  W,  NoverNo) 

C 

PRINT  * 

PRINT  *,  "NoverNo  (Exact)  = 

PRINT  *,  "Integrand*Disk  (Flat-top)  = 
PRINT  *,  "Cox-Munk  - 
C 

END 


" ,  NoverNo 
",  RqS 
",  RpG 
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